Take_Home Exercise 1: Geospatial Analytics for Social Good: Application of Spatial and Spatio-temporal Point Patterns Analysis to discover the geographical distribution of Armed Conflict in Myanmar

Author

Jiale SO

Published

August 22, 2024

Modified

September 22, 2024

1.0 Introduction

The study of armed conflicts in Myanmar has gained critical importance in understanding the geographical distribution and intensity of violence across different regions. Myanmar’s complex ethnic composition and ongoing civil strife make it a unique case for geospatial analysis. This project aims to apply spatial and spatio-temporal point pattern analysis methods to uncover the patterns of armed conflict between January 2021 and June 2024.

By leveraging conflict data from the Armed Conflict Location & Event Data (ACLED) and geospatial tools, we will focus on visualizing and interpreting conflict density through heat maps, Kernel Density Estimation (KDE), and advanced spatio-temporal analysis.

Our analysis will focus on four types of conflict events:

  1. Battles,
  2. Explosion/Remote violence,
  3. Strategic developments,
  4. Violence against civilians,

with particular attention paid to quarterly patterns in conflict occurrence.

2.0 Setting Up The Environment & Dataset

2.1 Installing the required Packages

Key Packages Used in the Project:

  1. sf: Handles simple features in R, allowing for spatial data manipulation and analysis. It is crucial for reading and managing geospatial data like shapefiles (e.g., Myanmar’s administrative boundaries).

  2. raster: Used for raster-based spatial data manipulation, especially for working with raster maps, such as Kernel Density Estimation (KDE) results.

  3. spatstat: A powerful package for spatial point pattern analysis. It helps to analyze and visualize spatial point data, particularly for identifying clusters or patterns in armed conflict events.

  4. sparr: Builds on spatstat and focuses on performing spatial and spatio-temporal kernel smoothing, which will be crucial for KDE and heatmap creation.

  5. tmap: A thematic mapping package that will allow us to create maps, including KDE visualizations on an OpenStreetMap base.

  6. tidyverse: A collection of data manipulation packages like dplyr, ggplot2, and purrr. It’s essential for data cleaning, manipulation, and visualization tasks.

  7. stpp: Used for spatio-temporal point pattern analysis, crucial for analyzing how conflict events evolve in both space and time.

  8. skimr: A quick and comprehensive tool to provide summaries and descriptive statistics for datasets, helping in the initial exploration.

  9. gganimate: Extends ggplot2 to create animated visualizations. We can use this for animated time-series or evolving conflict maps.

  10. ggplot2: The core plotting package in R, essential for creating visualizations like time series plots and KDE heatmaps.

  11. plotly: Useful for creating interactive visualizations, allowing users to explore spatial data interactively (e.g., hover over points to see conflict details).

  12. pacman: is a package management tool in R designed to streamline the process of loading and installing packages.

pacman::p_load(sf, raster, spatstat, sparr, tmap, tidyverse, stpp, skimr, gganimate, ggplot2, plotly, flexdashboard, DT,gridExtra, rlist, grid, animation)

2.2 Data-set involved in this topic

For this analysis, we use two key datasets:

2.2.1 ACLED Armed Conflict Data

Location & Event Data (ACLED)platform, which maintains an extensive record of conflict events globally. For this specific analysis, we will limit the dataset by filtering based on the following parameters to streamline data preparation and minimize the need for extensive data cleaning:

Data Parameter Filter Category
Date Range From 01/01/2021 to 30/06/2024.
Event Type 1. Battles
2. Violence Against Civilians
3. Explosions/Remote Violence
4. Strategic Developments
Country Myanmar
ACLED Configuration Image

Code to Import ACLED Dataset 
ACLEDData <- read_csv("data/raw/aspatial/2021-01-01-2024-06-30-Myanmar.csv")
Rows: 42608 Columns: 35
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (20): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (15): year, time_precision, inter1, inter2, interaction, iso, latitude, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

2.2.1.1 Understanding the data set fields.

Referencing this ACLED Official codebook, this is the dataset that we are working with, not to bore you with the details are mainly interested in the following fields,

  • Event ID: Unique identifier for each conflict event.

  • Event Date: Date of occurrence.

  • Event Type: Type of conflict event (e.g., Battles, Remote Violence).

  • Latitude/Longitude: Coordinates of the event.

  • Fatalities: Number of fatalities resulting from the event.

  • Actors: The groups or individuals involved in the conflict (e.g., state actors, ethnic armed groups).

  • Admin Levels: Administrative region, district, and township where the event took place.

If you’re interested in the data set fields to explore more, here’s the full fields!

ACLED Full Table Fields Summary
Fields name Fields Description Values
event_id_cnty A unique alphanumeric event identifier by number and country acronym. This identifier remains constant even when the event details are updated. E.g. ETH9766
event_date The date on which the event took place. Recorded as Year-Month-Day. E.g. 2023-02-16
year The year in which the event took place. E.g. 2018
time_precision A numeric code between 1 and 3 indicating the level of precision of the date recorded for the event. The higher the number, the lower the precision. 1, 2, or 3; with 1 being the most precise.
disorder_type The disorder category an event belongs to. Political violence, Demonstrations, or Strategic developments.
event_type The type of event; further specifies the nature of the event. E.g. BattlesFor the full list of ACLED event types, see the ACLED Event Types table.
sub_event_type A subcategory of the event type. E.g. Armed clashFor the full list of ACLED sub-event types, see the ACLED Event Types table.
actor1 One of two main actors involved in the event (does not necessarily indicate the aggressor). E.g. Rioters (Papua New Guinea)
assoc_actor_1 Actor(s) involved in the event alongside ‘Actor 1’ or actor designations that further identify ‘Actor 1’. E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank.
inter1 A numeric code between 0 and 8 indicating the type of ‘Actor 1’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). 1, 2, 3, 4, 5, 6, 7, or 8.
actor2 One of two main actors involved in the event (does not necessarily indicate the target or victim). E.g. Civilians (Kenya)Can be blank.
assoc_actor_2 Actor(s) involved in the event alongside ‘Actor 2’ or actor designation further identifying ‘Actor 2’. E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank.
inter2 A numeric code between 0 to 8 indicating the type of ‘Actor 2’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). 0, 1, 2, 3, 4, 5, 6, 7, or 8.
interaction A two-digit numeric code (combination of ‘Inter 1’ and ‘Inter 2’) indicating the two actor types interacting in the event (for more, see the section Actor Names, Types, and ‘Inter’ Codes). E.g.3, 58
civilian_targeting This column indicates whether the event involved civilian targeting. Either ‘Civilians targeted’ or blank.
iso A unique three-digit numeric code assigned to each country or territory according to ISO 3166. E.g. 231 for Ethiopia
region The region of the world where the event took place. E.g. Eastern Africa
country The country or territory in which the event took place. E.g. Ethiopia
admin1 The largest sub-national administrative region in which the event took place. E.g. Oromia
admin2 The second largest sub-national administrative region in which the event took place. E.g. ArsiCan be blank.
admin3 The third largest sub-national administrative region in which the event took place. E.g. MertiCan be blank.
location The name of the location at which the event took place. E.g. Abomsa
latitude The latitude of the location in four decimal degrees notation (EPSG:32647). E.g. 8.5907
longitude The longitude of the location in four decimal degrees notation (EPSG:32647). E.g. 39.8588
geo_precision A numeric code between 1 and 3 indicating the level of certainty of the location recorded for the event. The higher the number, the lower the precision. 1, 2, or 3; with 1 being the most precise.
source The sources used to record the event. Separated by a semicolon. E.g. Ansar Allah; Yemen Data Project
source_ scale An indication of the geographic closeness of the used sources to the event (for more, see the section Source Scale). E.g. Local partner-National
notes A short description of the event. E.g. On 16 February 2023, OLF-Shane abducted an unidentified number of civilians after stopping a vehicle in an area near Abomsa (Merti, Arsi, Oromia). The abductees were traveling from Adama to Abomsa, Arsi.
fatalities The number of reported fatalities arising from an event. When there are conflicting reports, the most conservative estimate is recorded. E.g. 3No information on fatalities is recorded as 0 reported fatalities.
tags Additional structured information about the event. Separated by a semicolon. E.g. women targeted: politicians; sexual violence
timestamp An automatically generated Unix timestamp that represents the exact date and time an event was uploaded to the ACLED API. E.g. 1676909320

2.2.2 Myanmar Administrative Boundaries (Shapefiles):

Obtained through Geonode Mimu, this shapefile helps us to build the map and set the boundary zone of each district of myanmar. This dataset provides the geographical boundaries of Myanmar’s administrative divisions, from the national level down to the township level. It is essential for mapping conflict events to specific regions.

Code to Import Shapefile Dataset 
M_State_Sf <- st_read(dsn="data/raw/geospatial/stateLevel", layer = "mmr_polbnda_adm1_250k_mimu_1") 
Reading layer `mmr_polbnda_adm1_250k_mimu_1' from data source 
  `C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial\stateLevel' 
  using driver `ESRI Shapefile'
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
M_State_Sf
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID          ST ST_PCODE           ST_RG          ST_MMR PCode_V
1         1  Ayeyarwady   MMR017          Region  ဧရာဝတီတိုင်းဒေသကြီး     9.4
2         2        Bago   MMR111          Region    ပဲခူးတိုင်းဒေသကြီး     9.4
3         4        Chin   MMR004           State       ချင်းပြည်နယ်     9.4
4         5      Kachin   MMR001           State       ကချင်ပြည်နယ်     9.4
5         6       Kayah   MMR002           State       ကယားပြည်နယ်     9.4
6         7       Kayin   MMR003           State        ကရင်ပြည်နယ်     9.4
7         8      Magway   MMR009          Region   မကွေးတိုင်းဒေသကြီး     9.4
8         9    Mandalay   MMR010          Region မန္တလေးတိုင်းဒေသကြီး     9.4
9        10         Mon   MMR011           State         မွန်ပြည်နယ်     9.4
10       11 Nay Pyi Taw   MMR018 Union Territory        နေပြည်တော်     9.4
                         geometry
1  MULTIPOLYGON (((95.20798 15...
2  MULTIPOLYGON (((96.17964 19...
3  MULTIPOLYGON (((93.36931 24...
4  MULTIPOLYGON (((97.59674 28...
5  MULTIPOLYGON (((97.1759 19....
6  MULTIPOLYGON (((97.81508 16...
7  MULTIPOLYGON (((94.11699 22...
8  MULTIPOLYGON (((96.14023 23...
9  MULTIPOLYGON (((97.73689 15...
10 MULTIPOLYGON (((96.32013 20...
M_District_Sf <- st_read(dsn="data/raw/geospatial", layer = "mmr_polbnda_adm2_250k_mimu") 
Reading layer `mmr_polbnda_adm2_250k_mimu' from data source 
  `C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
M_District_Sf
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID          ST ST_PCODE         DT   DT_PCODE      DT_MMR PCode_V
1         1  Ayeyarwady   MMR017   Hinthada MMR017D002    ဟင်္သာတခရိုင်     9.4
2         2  Ayeyarwady   MMR017    Labutta MMR017D004    လပွတ္တာခရိုင်     9.4
3         3  Ayeyarwady   MMR017     Maubin MMR017D005     မအူပင်ခရိုင်     9.4
4         4  Ayeyarwady   MMR017  Myaungmya MMR017D003 မြောင်းမြခရိုင်     9.4
5         5  Ayeyarwady   MMR017    Pathein MMR017D001      ပုသိမ်ခရိုင်     9.4
6         6  Ayeyarwady   MMR017     Pyapon MMR017D006     ဖျာပုံခရိုင်     9.4
7         7 Bago (East)   MMR007       Bago MMR007D001      ပဲခူးခရိုင်     9.4
8         8 Bago (East)   MMR007    Taungoo MMR007D002    တောင်ငူခရိုင်     9.4
9         9 Bago (West)   MMR008       Pyay MMR008D001      ပြည်ခရိုင်     9.4
10       10 Bago (West)   MMR008 Thayarwady MMR008D002   သာယာဝတီခရိုင်     9.4
                         geometry
1  MULTIPOLYGON (((95.12637 18...
2  MULTIPOLYGON (((95.04462 15...
3  MULTIPOLYGON (((95.38231 17...
4  MULTIPOLYGON (((94.6942 16....
5  MULTIPOLYGON (((94.27572 15...
6  MULTIPOLYGON (((95.20798 15...
7  MULTIPOLYGON (((95.90674 18...
8  MULTIPOLYGON (((96.17964 19...
9  MULTIPOLYGON (((95.70458 19...
10 MULTIPOLYGON (((95.85173 18...

2.2.2.1 Understanding the data set fields.

Field Name Description
OBJECTID This is a unique identifier for each feature in the dataset, typically used to identify individual records or polygons in the shapefile.
ST This represents the State or Region in Myanmar. For example, in your dataset, “Ayeyarwady” refers to a state/region.
ST_PCODE This stands for State Postal Code. It is a standardized code that represents each state or region in Myanmar, such as “MMR017” for Ayeyarwady.
DT This stands for District or Township within the respective state/region. For example, “Hinthada” is a district or township within Ayeyarwady.
DT_PCODE This stands for District/Township Postal Code. It is a standardized postal code for each district or township, such as “MMR017D002” for the Hinthada district/township in Ayeyarwady.
DT_MMR This field could be the District/Township name in Myanmar script, written in the local language. It may be an alternative representation of the “DT” field, showing the name of the district/township in Myanmar’s native language.
PCODE_V This could be a Version of the Postal Code or a verification value used internally in the dataset. In this case, the value is “9.4”, possibly indicating a specific version of postal codes or an accuracy measure.
geometry This column represents the spatial data for each district/township. It contains the geometrical shape (MULTIPOLYGON) defining the boundaries of the state or district/township, with coordinates provided in longitude and latitude.

3.0 Data Pre-processing

To ensure accuracy and usability of the data, several preprocessing steps will be undertaken for the different datasets.

3.1 Myanmar Shapefile

3.1.1 Setting the CRS for the

Since Myanmar uses CRS of 32647 and when we download the map it’s in WGS84, we should change it to 32647 .

st_crs(M_State_Sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_State_Sf <- st_transform(M_State_Sf, crs = 32647)
# Verify that the CRS has been correctly set
print(st_crs(M_State_Sf))
Coordinate Reference System:
  User input: EPSG:32647 
  wkt:
PROJCRS["WGS 84 / UTM zone 47N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 47N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",99,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
        BBOX[0,96,84,102]],
    ID["EPSG",32647]]
st_crs(M_District_Sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_District_Sf <- st_transform(M_District_Sf, crs = 32647)
# Verify that the CRS has been correctly set
print(st_crs(M_District_Sf))
Coordinate Reference System:
  User input: EPSG:32647 
  wkt:
PROJCRS["WGS 84 / UTM zone 47N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 47N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",99,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
        BBOX[0,96,84,102]],
    ID["EPSG",32647]]

3.1.2 Renaming and removal of column names

colnames(M_State_Sf) <- c("OBJECTID", "state","state_code","type", "state_myr", "mimi_version", "geometry")

M_State_Sf_Cleansed <- M_State_Sf %>% select(state, type, state_myr ,geometry)
colnames(M_District_Sf) <- c("OBJECTID", "state", "state_code", "district", "district_code", "district_mmr", "mimi_version", "geometry")

M_District_Sf_Cleansed <- M_District_Sf %>% select(state, district, district_mmr, geometry)

3.1.3 Checking for validity of data

When working with spatial data, it’s crucial to ensure that all geometries are valid. Invalid geometries can cause errors in analysis and visualization.

  1. Checking Validity with st_is_valid():

  2. Identifying Invalid Geometries:

  3. Fixing Invalid Geometries with st_make_valid()

#checking if it's valid
M_State_Sf_Validity <- st_is_valid(M_State_Sf_Cleansed)
M_State_Sf_Invalid <- which(!M_State_Sf_Validity)
if (length(M_State_Sf_Invalid) > 0) {
  print("MPZ Invalid!")
  print(M_State_Sf_cleansed[M_State_Sf_Invalid, ])
} else {
  print("it's valid!")
}
[1] "it's valid!"
#checking if it's valid
M_District_Sf_Validity <- st_is_valid(M_District_Sf_Cleansed)
M_District_Sf_Invalid <- which(!M_District_Sf_Validity)
if (length(M_District_Sf_Invalid) > 0) {
  print("MPZ Invalid!")
  print(M_District_Sf_cleansed[M_District_Sf_Invalid, ])
} else {
  print("it's valid!")
}
[1] "it's valid!"

3.1.4 Visualizing the mynamar map

On the top right, you can toggle between the district level and also the state level to understand more about the boundaries of Myanmar.

tm_shape(M_State_Sf_Cleansed) +  # Base map (Myanmar boundaries)
  tm_basemap(server = "OpenStreetMap.HOT") +  # Add OpenStreetMap as the basemap
  tm_polygons("state",  # Color the base map by state
              palette = "Set3",  # Use Set3 color palette for states
              border.col = "gray",  # Border color for the states
              alpha = 0.5,  # Semi-transparent polygons
              title = "State",  # Legend title for states
              legend.show = TRUE  # Show legend for state colors
             ) +
  tm_layout(main.title = "States of Myanmar",  # Main map title
            legend.outside = TRUE)  # Position the legend outside the map

There is more than 80 district here, so it only showcases 30 :)

tm_shape(M_District_Sf_Cleansed) +  # Base map (Myanmar boundaries)
  tm_basemap(server = "OpenStreetMap.HOT") +  # Add OpenStreetMap as the basemap
  tm_polygons("district",  # Color the base map by state
              palette = "Set3",  # Use Set3 color palette for states
              border.col = "gray",  # Border color for the states
              alpha = 0.5,  # Semi-transparent polygons
              title = "District",  # Legend title for states
              legend.show = TRUE  # Show legend for state colors
             ) +
  tm_layout(main.title = "District of Myanmar",  # Main map title
            legend.outside = TRUE)  # Position the legend outside the map
Warning: Number of levels of the variable "district" is 80, which is larger
than max.categories (which is 30), so levels are combined. Set
tmap_options(max.categories = 80) in the layer function to show all levels.
Some legend labels were too wide. These labels have been resized to 0.51, 0.48, 0.53, 0.52, 0.49. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

3.2 ACLED Data

3.2.1 Changing the Column Names

Since Myanmar’s regional hierarchy follows State, District, and Township levels, we will rename the columns accordingly:

  • admin1 → State

  • admin2 → District

  • admin3 → Township

This is important because different countries use different administrative hierarchies. For example, in Singapore, the hierarchy is organized by Region and Subzones. Adjusting these names ensures that our dataset aligns with Myanmar’s specific regional structure for accurate analysis.

ACLEDData_Cleanse <- ACLEDData %>%
  select(event_id_cnty, event_date, year, disorder_type, event_type, actor1, inter1, 
         actor2, inter2, interaction, admin1, admin2, admin3, location, latitude, 
         longitude, fatalities) %>%
  rename(state = admin1, district = admin2, township = admin3) %>%
  mutate(across(where(is.character), ~ replace_na(.x, "NA")),  # Replace NA in character columns with "NA"
         across(where(is.numeric), ~ replace_na(.x, 0)))  # Replace NA in numeric columns with 0

3.2.2 Adding a “Quarter-Year” Column

To facilitate our temporal analysis, we need to add a “quarter-year” column based on the event_date field. This can be done by adjusting the date format to represent the quarter and year, ensuring that each event is categorized by the specific quarter it occurred in (e.g., Q1-2021, Q2-2022). This will allow for easier analysis of conflict trends over time.

# Convert event_date to Date format (if it's not already a date)
ACLEDData_Cleanse$event_date <- as.Date(ACLEDData_Cleanse$event_date, format="%d-%b-%y")  # Adjust the format if needed
# Add a new column that shows the quarter and year
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  mutate(quarter_year = paste0("Q", quarter(event_date), "-", year(event_date)))

head(ACLEDData_Cleanse)
# A tibble: 6 × 18
  event_id_cnty event_date  year disorder_type   event_type actor1 inter1 actor2
  <chr>         <date>     <dbl> <chr>           <chr>      <chr>   <dbl> <chr> 
1 MMR64313      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
2 MMR64320      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
3 MMR64321      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
4 MMR64322      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
5 MMR64323      2024-06-30  2024 Political viol… Battles    PKDF …      3 Milit…
6 MMR64324      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
# ℹ 10 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
#   district <chr>, township <chr>, location <chr>, latitude <dbl>,
#   longitude <dbl>, fatalities <dbl>, quarter_year <chr>

3.2.3 Joining ACLED’s Codebook Description

ACLED’s stores their data for the column “interaction” and “inter1” and “inter2” in codes, using their code book, let’s reorganise their data for simplier view, we can reference the code book here to know what each code represent. Map it out as a csv file read it in and change accordingly.

3.2.3.1 Left joining inter1 and inter’s description.

For more details about each inter code read here.

ACLEDActorInterCode <- read_csv("data/raw/aspatial/ActorTypesInterCode.csv")
Rows: 8 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  left_join(ACLEDActorInterCode, by = c("inter1" = "code")) %>%
  rename(inter1_description = Description)
# Left join again for inter2
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  left_join(ACLEDActorInterCode, by = c("inter2" = "code")) %>%
  rename(inter2_description = Description)

head(ACLEDData_Cleanse)
# A tibble: 6 × 20
  event_id_cnty event_date  year disorder_type   event_type actor1 inter1 actor2
  <chr>         <date>     <dbl> <chr>           <chr>      <chr>   <dbl> <chr> 
1 MMR64313      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
2 MMR64320      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
3 MMR64321      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
4 MMR64322      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
5 MMR64323      2024-06-30  2024 Political viol… Battles    PKDF …      3 Milit…
6 MMR64324      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
# ℹ 12 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
#   district <chr>, township <chr>, location <chr>, latitude <dbl>,
#   longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
#   inter1_description <chr>, inter2_description <chr>

3.2.3.2 Left joining interaction description.

For more details about each interaction code read here.

ACLEDInteractionCode <- read_csv("data/raw/aspatial/AcledInteractionCodes.csv")
Rows: 44 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  left_join(ACLEDInteractionCode, by = c("interaction" = "code")) %>%
  rename(interaction_description = Description)
head(ACLEDData_Cleanse)
# A tibble: 6 × 21
  event_id_cnty event_date  year disorder_type   event_type actor1 inter1 actor2
  <chr>         <date>     <dbl> <chr>           <chr>      <chr>   <dbl> <chr> 
1 MMR64313      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
2 MMR64320      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
3 MMR64321      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
4 MMR64322      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
5 MMR64323      2024-06-30  2024 Political viol… Battles    PKDF …      3 Milit…
6 MMR64324      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
#   district <chr>, township <chr>, location <chr>, latitude <dbl>,
#   longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
#   inter1_description <chr>, inter2_description <chr>,
#   interaction_description <chr>

3.2.3 Making it a SF Object and reverse geolocate state and district

Since ACLED provides longitude and latitude data, I prefer to reverse geolocate the points to match Myanmar’s official state and district boundaries. We are uncertain how ACLED assigns these regions, so to ensure consistency, we remove the original state and district columns from the ACLED data and replace them with the geolocated values.

Steps:

  1. Convert ACLED Data to an SF Object: Using longitude and latitude coordinates, transform the ACLED dataset into a spatial object. Remember thatt we have to set CRS 32647 here as well.

  2. Perform a Spatial Join: Match the points from ACLED with the corresponding state and district boundaries from the m_sf shapefile, selecting only those columns.

  3. Remove Original Columns: After the spatial join, drop the original state, district, and township columns from the ACLED dataset.

  4. Rename the Joined Columns: Rename the newly joined state.y and district.y to state and district, effectively replacing the original columns with the reverse-geolocated values.

# Step 1: Convert ACLEDDataCleanse to an sf object using longitude and latitude
ACLEDData_Cleanse_Sf <- st_as_sf(ACLEDData_Cleanse, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

# Step 2: Transform ACLEDDataCleanse_sf to the same CRS as m_sf (EPSG: 32647)
ACLEDData_Cleanse_Sf <- st_transform(ACLEDData_Cleanse_Sf, crs = 32647)

# Step 3: Perform a spatial join, selecting only the state and district from m_sf
reverse_geolocated_state <- st_join(ACLEDData_Cleanse_Sf, M_State_Sf_Cleansed[, c("state")], join = st_intersects)

# Step 2: Spatial join to add 'district' from M_District_Sf_Cleansed
reverse_geolocated_district <- st_join(reverse_geolocated_state, M_District_Sf_Cleansed[, c("district")], join = st_intersects)

# Step 3: Remove original 'state', 'district', and 'township' columns from ACLEDData_Cleanse (if they exist)
# This step removes the original columns, and then renames the newly joined columns
ACLEDData_Cleanse_Sf <- reverse_geolocated_district %>%
  select(-contains("state.x"), -contains("district.x"), -contains("township")) %>%  # Remove old state, district, township columns
  rename(state = state.y, district = district.y)  # Rename newly joined columns

ACLEDData_Cleanse <- st_drop_geometry(ACLEDData_Cleanse_Sf)

# View the updated data
print(ACLEDData_Cleanse_Sf)
Simple feature collection with 42608 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -208804.4 ymin: 1103500 xmax: 640934.5 ymax: 3042960
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 42,608 × 21
   event_id_cnty event_date  year disorder_type  event_type actor1 inter1 actor2
   <chr>         <date>     <dbl> <chr>          <chr>      <chr>   <dbl> <chr> 
 1 MMR64313      2024-06-30  2024 Political vio… Battles    Peopl…      3 Milit…
 2 MMR64320      2024-06-30  2024 Political vio… Battles    Peopl…      3 Milit…
 3 MMR64321      2024-06-30  2024 Political vio… Battles    Peopl…      3 Milit…
 4 MMR64322      2024-06-30  2024 Strategic dev… Strategic… Milit…      1 NA    
 5 MMR64323      2024-06-30  2024 Political vio… Battles    PKDF …      3 Milit…
 6 MMR64324      2024-06-30  2024 Strategic dev… Strategic… Milit…      1 NA    
 7 MMR64325      2024-06-30  2024 Political vio… Battles    Milit…      1 PSLF/…
 8 MMR64326      2024-06-30  2024 Political vio… Battles    PSLF/…      2 Milit…
 9 MMR64328      2024-06-30  2024 Political vio… Battles    Milit…      1 PSLF/…
10 MMR64330      2024-06-30  2024 Political vio… Battles    Milit…      1 PSLF/…
# ℹ 42,598 more rows
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, location <chr>,
#   latitude <dbl>, longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
#   inter1_description <chr>, inter2_description <chr>,
#   interaction_description <chr>, geometry <POINT [m]>, state <chr>,
#   district <chr>

3.2.5 Visualizing it by Event Type

# Set tmap mode to plot for static maps
tmap_mode("view")
tmap mode set to interactive viewing
# Create the tmap object with the base map and event markers
tm_basemap(server = "OpenStreetMap") +  # Add OpenStreetMap as the base map
  
  # Add Myanmar state boundaries with transparency
  tm_shape(M_State_Sf_Cleansed) + 
  tm_polygons("state", alpha = 0.3, border.col = "gray", 
              title = "State Boundaries", legend.show = TRUE) +  # Add a legend for state boundaries
  
  # Add event markers (bubbles) with size based on fatalities
  tm_shape(ACLEDData_Cleanse_Sf) + 
  tm_bubbles(size = "fatalities",  # Marker size based on fatalities
             col = "event_type",  # Color markers by event type
             palette = "Set1",  # Use Set1 color palette for event types
             border.col = "black",  # Border color for bubbles
             border.alpha = 0.5,  # Semi-transparent border
             title.size = "Number of Fatalities",  # Title for bubble size legend
             title.col = "Event Types",  # Title for event type legend
             legend.size.show = TRUE,  # Show legend for bubble size
             legend.col.show = TRUE) +  # Show legend for event types
  
  # Layout settings for the map, including title and legend positioning
  tm_layout(main.title = "Myanmar's State Conflicts by Fatalities",  # Main map title
            main.title.size = 1.5,  # Title font size
            legend.outside = TRUE,  # Position legend outside the map
            legend.outside.size = 0.5,  # Adjust size of outside legend
            legend.position = c("left", "top"),  # Position for the event type legend
            legend.title.size = 1.2,  # Size of the legend title
            legend.text.size = 1,  # Size of the legend text
            legend.bg.color = "white",  # Background color for the legend
            legend.bg.alpha = 0.5,  # Transparency for the legend background
            inner.margins = c(0.05, 0.05, 0.05, 0.05))  # Inner margins for the map
legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
Legend for symbol sizes not available in view mode.

4.0 Exploratory Data analysis

With the data cleansed, let’s conduct a high-level analysis of our dataset to determine some key statistics.

The central question we’re exploring is:

“Which district in Myanmar has the highest number of conflicts that may be concerning to civilians?”

In section 3.2.5, we recognized that Myanmar has many states and numerous conflicts during this period. Therefore, we will aggregate the data at the state level to prioritize states based on the following key figures:

  • Civilian Involvement: We focus on events where either Actor 1 or Actor 2 is a civilian, as these are likely to raise humanitarian concerns.

  • Total Conflicts: This measures the overall number of conflict events in each state, helping us identify which states experienced the most conflict.

  • Conflict Density: This calculates the number of conflicts per square kilometer in each state, allowing us to understand the intensity of conflicts relative to the size of the state.

  • Total Fatalities: This sums up the total number of fatalities within each state, highlighting areas with the highest death toll, which are likely experiencing the most severe impacts.

  • Fatalities Density: This measures the fatalities per square kilometer, giving a sense of how deadly the conflicts are in each state in relation to its area.

This analysis will allow us to focus on the states most affected by conflicts that pose the greatest risk to civilians.

4.1 Aggregation of Data for exploratory purposes

Step 1: Filter Data for Civilians

In this first step, the data (ACLEDData_Cleanse) is filtered so that only rows where either “Civilians” are involved in the conflict (inter1_description or inter2_description) are kept. Additionally, rows where the state information is missing (NA) are removed because these conflicts can’t be assigned to a specific geographical area.

filtered_data <- ACLEDData_Cleanse %>%
  filter((inter1_description == "Civilians" | inter2_description == "Civilians") & !is.na(state))

Step 2: Calculate State Area in Square Kilometers

Next, the code calculates the area for each state from the M_State_Sf_Cleansed spatial dataset. The st_area function retrieves the area in square meters, and dividing by 1 million (1e6) converts it into square kilometers. This information will be used later to calculate conflict and fatality density.

state_area_km2 <- st_area(M_State_Sf_Cleansed) / 1e6
state_area_df <- data.frame(state = M_State_Sf_Cleansed$state, area_km2 = as.numeric(state_area_km2))

Step 3: Aggregate Total Conflicts and Fatalities by State

The next step is to group the filtered data by state and calculate two key metrics: the total number of conflicts and the total number of fatalities for each state. This is done using group_by and summarise.

agg_data_by_state <- filtered_data %>%
  group_by(state) %>%
  summarise(
    total_conflicts = n(), 
    total_fatalities = sum(fatalities, na.rm = TRUE), 
    .groups = "drop"
  )

Step 4: Convert to Non-Spatial Data Frame

Here, the spatial data is converted into a regular data frame by dropping the geometry information. This allows easier manipulation of the non-spatial data for further processing.

filtered_data_df <- st_drop_geometry(filtered_data)

Step 5: Prepare Data Frame for Yearly Data

To store yearly conflict and fatality data for each state, an empty data frame is initialized by selecting distinct states from the filtered data

final_result <- filtered_data_df %>% select(state) %>% distinct()

Step 6: Loop Through Each Year and Calculate Yearly Conflicts/Fatalities

In this step, the code iterates over each unique year in the dataset. For each year, it calculates the total conflicts and fatalities for each state. The results are merged back into the final_result data frame, with columns named according to the year (e.g., conflicts_2021, fatalities_2021).

unique_years <- unique(filtered_data_df$year)

for (yr in unique_years) {
  yearly_data <- filtered_data_df %>%
    filter(year == yr) %>%
    group_by(state) %>%
    summarise(
      conflicts = n(), 
      fatalities = sum(fatalities, na.rm = TRUE)
    ) %>%
    rename_at(vars(conflicts, fatalities), ~ paste0(., "_", yr))  # Rename with year
  
  final_result <- left_join(final_result, yearly_data, by = "state")
}

Step 7: Replace NAs with 0

Any missing data for a state-year combination (e.g., a state with no conflicts in a particular year) is replaced with 0 to avoid leaving gaps in the analysis.

final_result <- final_result %>%
  mutate(across(everything(), ~ ifelse(is.na(.), 0, .)))

Step 8: Calculate Conflict and Fatality Density

Now, the code calculates the density of conflicts and fatalities per square kilometer for each state. The results from Step 3 (agg_data_by_state) are joined with the state area data from Step 2. The density is simply the number of conflicts or fatalities divided by the area in square kilometers.

agg_data_with_density <- agg_data_by_state %>%
  left_join(state_area_df, by = "state") %>%
  mutate(
    conflict_density = total_conflicts / area_km2,
    fatality_density = total_fatalities / area_km2
  )

Step 9: Merge Yearly Data and Density with Spatial Data

Finally, the spatial data (M_State_Sf_Cleansed) is merged with the yearly conflict/fatality data (final_result) and the density data (agg_data_with_density). Additional calculations are made to find the fatality-per-conflict ratio which are then added to the final spatial dataset.

finalized_map <- M_State_Sf_Cleansed %>%
  left_join(final_result, by = "state") %>%
  left_join(agg_data_with_density, by = "state")

finalized_map <- finalized_map %>%
  mutate(conflict_fatality_ratio = ifelse(total_fatalities == 0, NA, total_fatalities / total_conflicts))

finalized_map <- finalized_map %>%
  mutate(
    fatality_per_conflict = ifelse(total_conflicts == 0, NA, total_fatalities / total_conflicts),
  )

finalized_map <- finalized_map %>%
  mutate(
    rank_total_conflicts = rank(-total_conflicts, ties.method = "min"),  # Rank total_conflicts (largest to smallest)
    rank_total_fatalities = rank(-total_fatalities, ties.method = "min"),  # Rank total_fatalities (largest to smallest)
    rank_conflict_density = rank(-conflict_density, ties.method = "min"),  # Rank conflict_density (largest to smallest)
    rank_fatality_density = rank(-fatality_density, ties.method = "min"),  # Rank fatality_density (largest to smallest)
    rank_fatality_per_conflict = rank(-fatality_per_conflict, ties.method = "min")  # Rank fatality_per_conflict (largest to smallest)
  )

Final Output: Display the Data

The finalized_map contains all the processed information, including conflict/fatality counts, densities, and state geometry, which can now be visualized or further analyzed.

head(finalized_map)
Simple feature collection with 6 features and 23 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -162806.6 ymin: 1682875 xmax: 491639.1 ymax: 3158467
Projected CRS: WGS 84 / UTM zone 47N
       state   type      state_myr conflicts_2024 fatalities_2024
1 Ayeyarwady Region ဧရာဝတီတိုင်းဒေသကြီး             64               9
2       Bago Region   ပဲခူးတိုင်းဒေသကြီး            194             128
3       Chin  State      ချင်းပြည်နယ်             72              52
4     Kachin  State      ကချင်ပြည်နယ်            160              95
5      Kayah  State      ကယားပြည်နယ်             21              38
6      Kayin  State       ကရင်ပြည်နယ်             52              59
  conflicts_2023 fatalities_2023 conflicts_2022 fatalities_2022 conflicts_2021
1             83               7            177              29            233
2            339             220            172              61            191
3            101              76            204              80            160
4            316             133            226             140            221
5            113              74            149              87            114
6            165              83            165              86            116
  fatalities_2021 total_conflicts total_fatalities area_km2 conflict_density
1              34             557               79 33738.82       0.01650917
2              87             896              496 38556.99       0.02323833
3              32             537              240 36276.56       0.01480295
4              71             923              439 88977.32       0.01037343
5             108             397              307 11763.10       0.03374961
6              57             498              285 30337.59       0.01641528
  fatality_density                       geometry conflict_fatality_ratio
1      0.002341516 MULTIPOLYGON (((93411.72 17...               0.1418312
2      0.012864076 MULTIPOLYGON (((203949.9 21...               0.5535714
3      0.006615843 MULTIPOLYGON (((-72918.03 2...               0.4469274
4      0.004933841 MULTIPOLYGON (((362696.3 31...               0.4756230
5      0.026098563 MULTIPOLYGON (((309155.7 22...               0.7732997
6      0.009394287 MULTIPOLYGON (((373551.3 18...               0.5722892
  fatality_per_conflict rank_total_conflicts rank_total_fatalities
1             0.1418312                   11                    14
2             0.5535714                    9                     6
3             0.4469274                   12                    13
4             0.4756230                    8                     7
5             0.7732997                   14                    10
6             0.5722892                   13                    12
  rank_conflict_density rank_fatality_density rank_fatality_per_conflict
1                    10                    15                         15
2                     9                     8                          4
3                    12                    11                          9
4                    15                    13                          6
5                     6                     4                          1
6                    11                    10                          3

4.2 Visualizing the aggregated data.

For each column that we created earlier, we can now showcase it in a choropleth map to highlight the states with the highest values in each category. These maps will help visually identify which state has the highest value for each column. You can navigate through the tabset to explore the different metrics and easily view the data on a map!

# Create the tmap object for total civilian conflicts
tm_total_conflicts <- 
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("total_conflicts", 
              title = "Total Civilian Related Conflicts",
              palette = "Reds", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar's State by Total Civilian Related Conflicts",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.2,  # Size of the legend title
    legend.text.size = 1,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )
tm_total_conflicts

# Create the tmap object for total civilian conflicts
tm_total_fatalities <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("total_fatalities", 
              title = "Total Fatalities by Civilian Related Conflicts",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar's State by Total Fatatlies For Civilian Related Conflicts",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.2,  # Size of the legend title
    legend.text.size = 1,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_total_fatalities

tm_conflict_density <-
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("conflict_density", 
              title = "Myanmar State by Civilian Related Conflicts Density",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar State by Civilian Related Conflicts Density",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.2,  # Size of the legend title
    legend.text.size = 1,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_conflict_density

tm_fatalities_conflict <- 
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("fatality_density", 
              title = "Myanmar State by Fatalies For Civilians Related Conflicts Density Per Km^2",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar State by Fatalies For Civilians Related Conflicts Density Per Km^2",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.0,  # Size of the legend title
    legend.text.size = 0.8,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_fatalities_conflict

tm_fatalities_density <-
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("fatality_per_conflict", 
              title = "Myanmar State by Fatalies for Civilian Related Conflicts",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar State by Fatalies for Civilian Related Conflicts",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.0,  # Size of the legend title
    legend.text.size = 0.8,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_fatalities_density

4.3 Analysis on the aggregated data to find top 3 state to analyse

4.3.1 Choose which States to analyse?

With the maps above, everyone can choose what they would like to analyze, whether it’s total civilian conflicts or fatalities by state. However, I would like to focus on the central question:

Which state has the highest amount of conflict and fatality density?

By using density, I’ve aimed to assess fatalities in relation to conflicts and ensure a fair comparison by adjusting for density.

Note

Keep in mind, this is not the most precise analysis, as other factors—such as ethnicity, population demographics, and gender—could also be considered. However, for the purpose of this study, I have chosen to focus on solely on the state’s conflict and fatality density.

Conflict Density: This is the number of conflicts per square kilometer in a state.

  1. Conflict Density=Total ConflictsArea of State (km²) = Conflict Density=Area of State (km²)Total Conflicts​

    This tells you how many conflicts occur per unit of area (1 km²) in each state.

  2. Fatality Density: This is the number of fatalities per square kilometer.

    Fatality Density=Total FatalitiesArea of State (km²) = Fatality Density=Area of State (km²)Total Fatalities​

4.3.2 Viewing the top 3 states across

Using the map below, you can interactively explore the Conflict and Fatality Density. We observe that the top three states are as follows, and here’s how they compare across the data set.

You can use the layer toggle to switch between fatalities and conflicts, or click on the map layers to view detailed, aggregated data for each state.

# Step 1: Set tmap mode to view for interactive maps

# Step 2: Define popup variables to show all relevant columns when clicking on a state
popup_variables <- c( 
                     "Total Civilian Conflicts" = "total_conflicts",
                     "Fatalities" = "total_fatalities",
                     "Civilian Conflict Density" = "conflict_density",
                     "Fatality Density" = "fatality_density",
                     "Fatality Per Civilian Conflict" = "fatality_per_conflict",
                     "Ranked For Total Conflict By State (Out of 15)" = "rank_total_conflicts",
                     "Ranked For Total Fatalities By State (Out of 15)" = "rank_total_fatalities",
                     "Ranked For Conflict Density By State (Out of 15)" = "rank_conflict_density",
                     "Ranked For Fatalities Density By State (Out of 15)" = "rank_fatality_density",
                     "Ranked For Fatalies/Conflict By State (Out of 15)" = "rank_fatality_per_conflict"
                     )

# Step 3: Create the map with 6 layers for toggling
tm <- tm_basemap(server = "OpenStreetMap") +
  tm_shape(finalized_map) +
  tm_polygons("conflict_density", 
              title = "Civilian Conflict Density",
              palette = "Oranges", 
              border.col = "black",
              popup.vars = popup_variables,
              id = "state", 
              group = "Civilian Conflict Density") +
  tm_shape(finalized_map) +
  tm_polygons("fatality_density", 
              title = "Fatality Density",
              palette = "Purples", 
              border.col = "black",
              popup.vars = popup_variables,
              id = "state", 
              group = "Fatality Density") +
  tm_layout(
    legend.outside = TRUE,
    legend.outside.size = 0.5,  # Adjust the size of the legend
    legend.position = c("left", "top")  # Position the legend
  )
# Step 5: Display the interactive map
tm
legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.

Here we build the bar chart to view rank the states base on the variables.

finalized_data_df <- as.data.frame(st_drop_geometry(finalized_map))
# Sort the data by conflict_density and fatality_density in descending order
finalized_data_df <- finalized_data_df %>%
  arrange(desc(conflict_density), desc(fatality_density))

# Format the densities to show "per km²"
finalized_data_df$conflict_density_label <- paste0(round(finalized_data_df$conflict_density, 2), " per km²")
finalized_data_df$fatality_density_label <- paste0(round(finalized_data_df$fatality_density, 2), " per km²")

# Create bar chart for Conflict Density (sorted by conflict_density)
conflict_density_plot <- ggplot(finalized_data_df, aes(x = reorder(state, -conflict_density), y = conflict_density)) +
  geom_bar(stat = "identity", fill = "orange") +
  geom_text(aes(label = conflict_density_label), vjust = 1.8, size = 2.0, color = "black", position = position_stack(vjust = 0.5)) +  # Place text at bottom of bar
  ggtitle("Top States by Conflict Density (per km²)") +
  xlab("State") +
  ylab("Conflict Density (per km²)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Create bar chart for Fatality Density (sorted by fatality_density)
fatality_density_plot <- ggplot(finalized_data_df, aes(x = reorder(state, -fatality_density), y = fatality_density)) +
  geom_bar(stat = "identity", fill = "purple") +
  geom_text(aes(label = fatality_density_label), vjust = 1.8, size = 2.0, color = "black", position = position_stack(vjust = 0.5)) +  # Place text at bottom of bar
  ggtitle("Top States by Fatality Density (per km²)") +
  xlab("State") +
  ylab("Fatality Density (per km²)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Display both charts side by side
grid.arrange(conflict_density_plot, fatality_density_plot, ncol = 2)

4.3.3 Top 3 Chosen state for further analysis

4.3.3.1 Yangon

Yangon, Myanmar’s largest city and former capital, remains the country’s economic and cultural hub. Although it is no longer the political capital, Yangon holds the highest population density in the country, making it a vital urban center with significant influence.

According to the data, Yangon ranks 5th for total conflicts, with 1,489 civilian conflicts reported, and 8th in fatalities, with 431 deaths. However, it stands out with the highest civilian conflict density (0.0531) and fatality density (0.01536) among Myanmar’s 15 states. This means that while the overall number of conflicts may not be the highest, they are more concentrated within Yangon’s urban area, resulting in the highest per-unit area conflict rate.

Yangon’s urbanization and strategic economic importance likely make it a hotspot for civilian unrest. Social tensions in a densely populated city lead to frequent conflicts, reflected in its rank of 1st for conflict density and fatality density. Despite fewer total conflicts compared to other states, Yangon’s density of unrest highlights its vulnerability to civil tensions.

4.3.3.2 Mandalay

Mandalay, the second-largest city in Myanmar, is known for its cultural significance and central location in the country. It holds a key strategic position in the economy and transportation networks, making it vulnerable to unrest.

The data shows Mandalay ranks 2nd for total civilian conflicts with 2,021 conflicts and 950 fatalities. Its civilian conflict density is 0.0255, placing it 3rd in density, while its fatality density is 0.0120, ranked 2nd. Despite being second in conflict counts, the concentration of conflicts and fatalities suggests that civil unrest is deeply rooted in the region, but not as dense as in Yangon.

Mandalay’s high fatality-to-conflict ratio (0.4701 fatalities per conflict) highlights the deadly nature of conflicts, ranking 7th in fatalities per conflict. This reflects a region where conflicts, though frequent, lead to significant loss of life compared to other states.

4.3.3.3 Sagaling.

Sagaing, located in northern Myanmar, is one of the country’s most conflict-ridden regions, often affected by armed insurgencies and ethnic clashes due to its proximity to volatile borders.

The data reveals that Sagaing ranks 1st in total civilian conflicts with 5,346 conflicts and 3,319 fatalities, giving it the highest fatality count as well. However, its civilian conflict density is 0.0241 (ranked 3rd) and its fatality density is 0.01497 (ranked 3rd), indicating a widespread yet slightly less concentrated conflict zone compared to Yangon and Mandalay.

Sagaing’s fatality per civilian conflict ratio of 0.6208 ranks 2nd, emphasizing the high lethality of conflicts in the region. With both the largest number of conflicts and a high fatality rate, Sagaing represents one of the most dangerous and unstable areas in Myanmar for civilians/

4.4 Final Round of Cleansing ACLED Data

In this final stage of the data cleansing process, we focus on preparing the ACLED (Armed Conflict Location & Event Data) for the different states. The key steps include:

  1. Filtering: We first filter the dataset to focus on data specific to Yangon.

  2. Converting to Spatial Format: Using the st_as_sf() function, we transform the filtered data into a simple features (sf) object, a format commonly used for spatial data in R. This involves specifying longitude and latitude as the coordinates and setting the correct WGS84 (CRS 4326) coordinate system for global geographic data.

  3. Coordinate System Transformation: We apply the st_transform() function to project the spatial data from WGS84 (CRS 4326) to UTM Zone 47N (CRS 32647), which is more appropriate for spatial analysis in Myanmar. This step ensures accurate distance and area calculations.

  4. Validity Check: Using st_is_valid(), we validate the geometries in the spatial data to ensure there are no invalid shapes, such as self-intersecting polygons, which could cause errors during analysis.

  5. Conversion to Spatial Data Frame: Finally, we convert the sf object into a Spatial Data Frame (as_Spatial()) for compatibility with specific spatial analysis functions that require this format.

This process ensures that the States data is accurately represented in a geospatial format, with proper coordinate transformation and validity checks applied. It prepares the data for further geospatial operations, such as mapping and spatial pattern analysis, with the confidence that the dataset is clean and valid for use.

# Filter data for Yangon state from the entire ACLED dataset
Yangon_ACLED_Data <- filtered_data %>% 
  filter(state %in% c("Yangon"))

# Check the class of the data to confirm it's a dataframe
class(Yangon_ACLED_Data)
[1] "tbl_df"     "tbl"        "data.frame"
# Convert the filtered data into a spatial format (simple features object)
# Use longitude and latitude columns as coordinates, with WGS84 (CRS 4326) as the coordinate system
Yangon_ACLED_Data_Sf <- st_as_sf(Yangon_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

# Transform the coordinate system to UTM Zone 47N (CRS 32647) for better spatial precision in Myanmar
Yangon_ACLED_Data_Sf <- st_transform(Yangon_ACLED_Data_Sf, crs = 32647)

# Check the class of the spatial object to confirm conversion
class(Yangon_ACLED_Data_Sf)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
# Validate the spatial object to ensure all geometries are valid (no broken or self-intersecting geometries)
Yangon_ACLED_Validity <- st_is_valid(Yangon_ACLED_Data_Sf)

# Identify any invalid geometries and print them if they exist
Yangon_invalid <- which(!Yangon_ACLED_Validity)
if (length(Yangon_invalid) > 0) {
  print("Yangon is Invalid!")
  print(Yangon_ACLED_Data_Sf[Yangon_invalid, ])
} else {
  print("Yangon_ACLED_Data_Sf is valid!")
}
[1] "Yangon_ACLED_Data_Sf is valid!"
# Convert the sf object into a Spatial Data Frame (for compatibility with certain spatial analysis functions)
Yangon_ACLED_SFDF <- as_Spatial(Yangon_ACLED_Data_Sf)

# Final output - Yangon ACLED data in spatial data frame format, ready for analysis
Yangon_ACLED_SFDF
class       : SpatialPointsDataFrame 
features    : 1489 
extent      : 148595, 259745.5, 1819680, 1944325  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=47 +datum=WGS84 +units=m +no_defs 
variables   : 20
names       : event_id_cnty, event_date, year,          disorder_type,                 event_type,                       actor1, inter1,                actor2, inter2, interaction,      location, latitude, longitude, fatalities, quarter_year, ... 
min values  :      MMR10949,      18633, 2021,     Political violence, Explosions/Remote violence,           5 Brothers Younger,      1, Civilians (Australia),      0,          17,     Ah Hpyauk,  16.4389,   95.6934,          0,      Q1-2021, ... 
max values  :      MMR64089,      19895, 2024, Strategic developments, Violence against civilians, YUG: Yangon Urban Guerrillas,      7,                    NA,      7,          70, Zee Hpyu Kone,  17.5612,   96.7448,         13,      Q4-2023, ... 
# Filter data for Mandalay state from the entire ACLED dataset
Mandalay_ACLED_Data <- filtered_data %>% 
  filter(state %in% c("Mandalay"))
# Check the class of the filtered data to confirm it's still a dataframe
class(Mandalay_ACLED_Data)
[1] "tbl_df"     "tbl"        "data.frame"
# Convert the filtered data into a spatial format (simple features object)
# Longitude and latitude are used as coordinates, and WGS84 (CRS 4326) is set as the coordinate system
Mandalay_ACLED_Data_Sf <- st_as_sf(Mandalay_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

# Transform the coordinate system to UTM Zone 47N (CRS 32647) for better accuracy in the Myanmar region
Mandalay_ACLED_Data_Sf <- st_transform(Mandalay_ACLED_Data_Sf, crs = 32647)

# Verify that the data has been successfully converted into a spatial object
class(Mandalay_ACLED_Data_Sf)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
# Validate the spatial data to check for any invalid geometries (e.g., self-intersections, broken polygons)
Mandalay_ACLED_Validity <- st_is_valid(Mandalay_ACLED_Data_Sf)

# Identify any invalid geometries and print them for inspection, if found
Mandalay_invalid <- which(!Mandalay_ACLED_Validity)
if (length(Mandalay_invalid) > 0) {
  print("Mandalay is Invalid!")
  print(Mandalay_ACLED_Data_Sf[Mandalay_invalid, ])
} else {
  # If all geometries are valid, print confirmation
  print("Mandalay_ACLED_Data_Sf is valid!")
}
[1] "Mandalay_ACLED_Data_Sf is valid!"
# Convert the simple features object to a Spatial Data Frame for compatibility with spatial functions
Mandalay_ACLED_SFDF <- as_Spatial(Mandalay_ACLED_Data_Sf)

# Final output - Mandalay ACLED data in spatial data frame format, ready for analysis
Mandalay_ACLED_SFDF
class       : SpatialPointsDataFrame 
features    : 2021 
extent      : 68129.21, 255317.2, 2261881, 2620357  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=47 +datum=WGS84 +units=m +no_defs 
variables   : 20
names       : event_id_cnty, event_date, year,          disorder_type,                 event_type,                                   actor1, inter1,              actor2, inter2, interaction,         location, latitude, longitude, fatalities, quarter_year, ... 
min values  :      MMR10960,      18640, 2021,     Political violence, Explosions/Remote violence, 21KPG: 21 Guerrilla Force - Kyaukpadaung,      1, Civilians (Myanmar),      0,          17, 4 Maing Kan Thar,  20.4319,    94.849,          0,      Q1-2021, ... 
max values  :      MMR64318,      19904, 2024, Strategic developments, Violence against civilians,          Zero Guerrilla Force - Myingyan,      7,                  NA,      7,          70,       Zin Chaung,   23.667,   96.6148,         10,      Q4-2023, ... 
Sagaing_ACLED_Data <- filtered_data %>% 
  filter(state %in% c("Sagaing"))

# Check the class of the filtered data to confirm it's still a dataframe
class(Sagaing_ACLED_Data)
[1] "tbl_df"     "tbl"        "data.frame"
# Convert the filtered data into a spatial format (simple features object)
# Longitude and latitude are used as coordinates, and WGS84 (CRS 4326) is set as the coordinate system
Sagaing_ACLED_Data_Sf <- st_as_sf(Sagaing_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

# Transform the coordinate system to UTM Zone 47N (CRS 32647) for better accuracy in the Myanmar region
Sagaing_ACLED_Data_Sf <- st_transform(Sagaing_ACLED_Data_Sf, crs = 32647)

# Verify that the data has been successfully converted into a spatial object
class(Sagaing_ACLED_Data_Sf)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
# Validate the spatial data to check for any invalid geometries (e.g., self-intersections, broken polygons)
Sagaing_ACLED_Validity <- st_is_valid(Sagaing_ACLED_Data_Sf)

# Identify any invalid geometries and print them for inspection, if found
Sagaing_invalid <- which(!Sagaing_ACLED_Validity)
if (length(Sagaing_invalid) > 0) {
  print("Sagaing is Invalid!")
  print(Sagaing_ACLED_Data_Sf[Sagaing_invalid, ])
} else {
  # If all geometries are valid, print confirmation
  print("Sagaing_ACLED_Data_Sf is valid!")
}
[1] "Sagaing_ACLED_Data_Sf is valid!"
# Convert the simple features object to a Spatial Data Frame for compatibility with spatial functions
Sagaing_ACLED_SFDF <- as_Spatial(Sagaing_ACLED_Data_Sf)

# Final output - Sagaing ACLED data in spatial data frame format, ready for analysis
Sagaing_ACLED_SFDF
class       : SpatialPointsDataFrame 
features    : 5346 
extent      : -16397.05, 256586.9, 2393568, 2914063  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=47 +datum=WGS84 +units=m +no_defs 
variables   : 20
names       : event_id_cnty, event_date, year,          disorder_type,                 event_type,                                      actor1, inter1,            actor2, inter2, interaction,               location, latitude, longitude, fatalities, quarter_year, ... 
min values  :      MMR10952,      18640, 2021,     Political violence, Explosions/Remote violence, ABSDF: All Burma Students' Democratic Front,      1, Civilians (China),      0,          17,               55 Maing,   21.605,   93.9575,          0,      Q1-2021, ... 
max values  :      MMR65891,      19898, 2024, Strategic developments, Violence against civilians,             Zero Guerrilla Force - Myingyan,      7,                NA,      7,          70, Zin Ka Le (Zin Ka Lin),  26.3006,   96.6034,        175,      Q4-2023, ...

5.0 Spatial Data Wrangling

Spatial data wrangling is a crucial step in preparing the dataset for spatial point pattern analysis. This section involves setting up the spatial boundaries for each state and transforming the data into appropriate formats for analysis. The goal is to define the area within which the point patterns (conflict events) will be analyzed and ensure that the data is in the correct format for further spatial operations.

5.1 Setting up the Owin Window for States.

In this step, we are setting up an owin object for each state using the spatstat package. An owin object defines the spatial window—or the boundary—within which our point pattern analysis will take place. Essentially, this window will encapsulate the polygonal boundary of each state, allowing us to focus on conflict data within those precise boundaries.

To create the owin object, we will:

  • Extract the polygon boundaries for each state from the spatial data.

  • Convert these boundaries into an owin object using the as.owin() function.

  • These spatial windows will ensure that all spatial operations, such as density estimation and pattern analysis, occur strictly within the boundaries of the specified state (e.g., Yangon, Mandalay, Sagaing).

This setup is crucial for ensuring that our spatial analysis focuses only on relevant geographic areas and doesn’t incorporate any points outside the boundaries.

Yangon_Sf <- M_State_Sf %>% 
  filter(state %in% c("Yangon"))
# Step 1: Extract the individual polygons from the multipolygon
yangon_polygons <- st_cast(Yangon_Sf$geometry, "POLYGON")
yangon_polygons_filtered <- yangon_polygons[-c(1, 2)]
yangon_multipolygon_filtered <- st_combine(yangon_polygons_filtered)
# Add the filtered multipolygon back to the Yangon_Sf object
Yangon_Sf$geometry <- yangon_multipolygon_filtered
Yangon <- as_Spatial(Yangon_Sf)
Yangon_Owin <- as.owin(Yangon_Sf)
plot(Yangon_Owin)

summary(Yangon_Owin)
Window: polygonal boundary
6 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1            3832  9.84832e+09      1.00e+00
polygon 2 (hole)        3 -1.78834e+02     -1.82e-08
polygon 3              14  2.11571e+05      2.15e-05
polygon 4              11  1.59536e+05      1.62e-05
polygon 5              20  1.79201e+06      1.82e-04
polygon 6              37  2.56010e+06      2.60e-04
enclosing rectangle: [140506.53, 268404.01] x [1805599.6, 1970174] units
                     (127900 x 164600 units)
Window area = 9853040000 square units
Fraction of frame area: 0.468
# Filter the dataset to only include data for Mandalay state
Mandalay_Sf <- M_State_Sf %>% 
  filter(state %in% c("Mandalay"))
# Convert the Mandalay spatial dataframe into a Spatial Data Frame for further spatial operations
Mandalay <- as_Spatial(Mandalay_Sf)

# Convert the Mandalay spatial dataframe into an owin object for spatial point pattern analysis
Mandalay_Owin <- as.owin(Mandalay_Sf)

# Plot the owin object to visualize the boundaries of Mandalay
plot(Mandalay_Owin)

# Provide a summary of the owin object, showing its properties like dimensions and bounding box
summary(Mandalay_Owin)
Window: polygonal boundary
single connected closed polygon with 5914 vertices
enclosing rectangle: [66291.97, 282396.46] x [2234807.4, 2622191.5] units
                     (216100 x 387400 units)
Window area = 30998600000 square units
Fraction of frame area: 0.37
# Filter the dataset to only include data for Sagaing state
Sagaing_Sf <- M_State_Sf %>% 
  filter(state %in% c("Sagaing"))

# Convert the Sagaing spatial dataframe into a Spatial Data Frame for further spatial operations
Sagaing <- as_Spatial(Sagaing_Sf)

# Convert the Sagaing spatial dataframe into an owin object for spatial point pattern analysis
Sagaing_Owin <- as.owin(Sagaing_Sf)
# Plot the owin object to visualize the boundaries of Sagaing
plot(Sagaing_Owin)

# Provide a summary of the owin object, showing its properties like dimensions and bounding box
summary(Sagaing_Owin)
Window: polygonal boundary
single connected closed polygon with 5882 vertices
enclosing rectangle: [-17699.96, 308341.37] x [2390344.6, 3029739.1] units
                     (326000 x 639400 units)
Window area = 9.3875e+10 square units
Fraction of frame area: 0.45

5.2 Setting up the ACLED Spatial Class

After defining the owin window for each state, we need to convert the conflict data into a suitable spatial point pattern object (typically using spatstat’s ppp class). This class is essential for performing spatial point pattern analysis, as it links the geographic coordinates of conflict events to the state boundary.

Here’s what we’ll do:

  • We will first ensure that the conflict data is transformed into the appropriate spatial projection (UTM or other applicable projections).

  • The conflict data points will be linked to the owin window to define where these events occur within the state’s boundary.

  • The result will be a ppp (point pattern) object, which is necessary for conducting spatial operations like density estimation, K-function analysis, or G-function analysis.

This step transforms our dataset into a fully spatially aware format, ready for statistical analysis.

# Extract the coordinates (longitude, latitude) from the Yangon spatial object
Yangon_ACLED_coords <- st_coordinates(Yangon_ACLED_Data_Sf)

# Define the bounding box (xmin, xmax, ymin, ymax) for Yangon, which sets the spatial extent of the window
Yangon_ACLED_bbox <- st_bbox(Yangon_ACLED_Data_Sf)

# Create the spatial window (owin object) for Yangon using the bounding box ranges
Yangon_ACLED_window <- owin(xrange = Yangon_ACLED_bbox[c("xmin", "xmax")], yrange = Yangon_ACLED_bbox[c("ymin", "ymax")])

# Create a ppp (point pattern) object for Yangon using the extracted coordinates and the defined window
Yangon_ACLED_ppp <- ppp(x = Yangon_ACLED_coords[, 1], y = Yangon_ACLED_coords[, 2], window = Yangon_ACLED_window)
# Check the summary of the ppp object to verify the number of points, window dimensions, and other properties
summary(Yangon_ACLED_ppp)
Planar point pattern:  1489 points
Average intensity 1.074751e-07 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 10 decimal places

Window: rectangle = [148595.03, 259745.55] x [1819680.1, 1944325.3] units
                    (111200 x 124600 units)
Window area = 13854400000 square units
# Plot the ppp object to visually inspect the spatial distribution of conflict points in Yangon
plot(Yangon_ACLED_ppp)

# Extract the coordinates (longitude, latitude) from the Mandalay spatial object
Mandalay_ACLED_coords <- st_coordinates(Mandalay_ACLED_Data_Sf)

# Define the bounding box (xmin, xmax, ymin, ymax) for Mandalay, which sets the spatial extent of the window
Mandalay_ACLED_bbox <- st_bbox(Mandalay_ACLED_Data_Sf)

# Create the spatial window (owin object) for Mandalay using the bounding box ranges
Mandalay_ACLED_window <- owin(xrange = Mandalay_ACLED_bbox[c("xmin", "xmax")], yrange = Mandalay_ACLED_bbox[c("ymin", "ymax")])

# Create a ppp (point pattern) object for Mandalay using the extracted coordinates and the defined window
Mandalay_ACLED_ppp <- ppp(x = Mandalay_ACLED_coords[, 1], y = Mandalay_ACLED_coords[, 2], window = Mandalay_ACLED_window)

# Check the summary of the ppp object to verify the number of points, window dimensions, and other properties
summary(Mandalay_ACLED_ppp)
Planar point pattern:  2021 points
Average intensity 3.011815e-08 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 11 decimal places

Window: rectangle = [68129.21, 255317.22] x [2261880.7, 2620356.6] units
                    (187200 x 358500 units)
Window area = 67102400000 square units
# Plot the ppp object to visually inspect the spatial distribution of conflict points in Mandalay
plot(Mandalay_ACLED_ppp)

# Extract the coordinates (longitude, latitude) from the Sagaing spatial object
Sagaing_ACLED_coords <- st_coordinates(Sagaing_ACLED_Data_Sf)

# Define the bounding box (xmin, xmax, ymin, ymax) for Sagaing, which sets the spatial extent of the window
Sagaing_ACLED_bbox <- st_bbox(Sagaing_ACLED_Data_Sf)

# Create the spatial window (owin object) for Sagaing using the bounding box ranges
Sagaing_ACLED_window <- owin(xrange = Sagaing_ACLED_bbox[c("xmin", "xmax")], yrange = Sagaing_ACLED_bbox[c("ymin", "ymax")])

# Create a ppp (point pattern) object for Sagaing using the extracted coordinates and the defined window
Sagaing_ACLED_ppp <- ppp(x = Sagaing_ACLED_coords[, 1], y = Sagaing_ACLED_coords[, 2], window = Sagaing_ACLED_window)
# Check the summary of the ppp object to verify the number of points, window dimensions, and other properties
summary(Sagaing_ACLED_ppp)
Planar point pattern:  5346 points
Average intensity 3.762491e-08 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 13 decimal places

Window: rectangle = [-16397.05, 256586.9] x [2393568.1, 2914062.9] units
                    (273000 x 520500 units)
Window area = 1.42087e+11 square units
# Plot the ppp object to visually inspect the spatial distribution of conflict points in Sagaing
plot(Sagaing_ACLED_ppp)

Note

Why dont you jitter, delete or mark any conflicts here?

The reason I did not apply jittering, deletion, or marking of the conflicts at this stage is that all the conflicts are unique and occurred across different time zones. Since there is no overlap in terms of exact location or time, it makes no sense to modify the data at this point. We will consider jittering later if we find that events are clustered too closely together in space, but for now, each conflict is distinct, and no adjustments are necessary.

5.3 Combining them both

Finally, we will merge the owin object (state boundary) and the spatial point pattern data (conflict events) into a single object. This combination ensures that:

  • The spatial window constrains the analysis, focusing on the area within each state’s boundary.

  • All conflict events are accurately represented within the confines of the defined boundary.

This step allows for robust spatial point pattern analysis, ensuring that both the boundaries and events are properly accounted for. With this merged object, we can perform advanced geospatial techniques, such as analyzing the distribution of conflict events, identifying hotspots, and calculating density functions within the state.

Yangon_ACLED_ppp = Yangon_ACLED_ppp[Yangon_Owin]
plot(Yangon_ACLED_ppp)

Mandalay_ACLED_ppp = Mandalay_ACLED_ppp[Mandalay_Owin]
plot(Mandalay_ACLED_ppp)

Sagaing_ACLED_ppp = Sagaing_ACLED_ppp[Sagaing_Owin]
plot(Sagaing_ACLED_ppp)

6.0 Summary of Data sets Prepared

We have created several distinct datasets, each serving a specific role in our spatial analysis of conflict events in Myanmar. Here’s a summary of each dataset, including the type of data they represent and a brief description of their purpose:

  1. M_State_SF (sf object)

    • Type: Simple feature object (sf)

    • Description: This is the cleaned and formatted dataset containing conflict events from ACLED (Armed Conflict Location & Event Data). It includes attributes like location, event type, and time period. This dataset serves as the foundation for further spatial analysis.

  2. {State}_SF (sf object)

    • Type: Simple feature object (sf)

    • Description: This dataset represents the Myanmar boundary and is used as the base map for the analysis. It defines the spatial extent and provides a reference layer for plotting conflict events within the country.

  3. Finalized_map (aggregated data)

    • Type: Aggregated dataset

    • Description: This dataset contains aggregated data for different event types, quarters, or other relevant categories. It is primarily used for exploratory analysis to get an overview of trends and patterns before performing detailed spatial point pattern analysis.

  4. filtered_data (filtered ACLED data)

    • Type: Filtered dataset

    • Description: This dataset consists of ACLED data filtered specifically for civilian-related conflicts. It focuses on events where civilians are involved, either in inter 1 or inter 2, and is the primary dataset for analyzing civilian conflict events.

  5. {State}_ACLED_Data_Sf (sf object)

    • Type: Simple feature object (sf)

    • Description: These datasets (e.g., Yangon_ACLED_Data_Sf, Mandalay_ACLED_Data_Sf) are state-specific ACLED data stored as sf objects. They contain conflict data for individual states and will be used to analyze conflict patterns within each state’s boundary.

  6. {State}_ACLED_ppp (ppp object)

    • Type: Planar point pattern object (ppp)

    • Description: These are point pattern datasets created from the state-specific ACLED data (e.g., Yangon_ACLED_ppp, Mandalay_ACLED_ppp). They store the coordinates of conflict events and are used for in-depth spatial analyses, such as density estimation and clustering. If needed, jittering is applied to handle close events.

  7. {State}_Owin(owin object)

    • Type: Spatial window object (owin)

    • Description: These objects represent the spatial boundary of each state (e.g., Yangon_Owin, Mandalay_Owin). The owin window defines the geographic limits for the spatial point pattern analysis, ensuring that all conflict events are analyzed within the correct state boundary.

7.0 Deriving the Quarterly Kernel Density Estimation Layers

Kernel Density Estimation (KDE) is a key technique for spatial analysis, allowing us to estimate the probability density of events occurring in space. By applying KDE across different quarters for each state (Yangon, Mandalay, and Sagaing), we can visualize how the intensity of conflict events evolves over time. This will give us insights into hotspots and help track the spread or concentration of conflicts across Myanmar. Our analysis will focus on generating KDE layers for each quarter, comparing them, and identifying significant patterns.

7.1 Setting up a function to create quarterly for each state

We will create a function that takes the conflict data for each state and computes the quarterly KDE. This function will:

  • Filter the conflict events by quarter.

  • Apply the KDE to estimate the spatial density of events for that quarter.

  • Return a KDE layer that can be used for visualization and analysis.

Full Code 
process_quarter_conflicts <- function(region_sf, region_window, region_name, data_sf, sigma_type, kernel_method) {
  # Extract unique quarters (reversed for correct order)
  quarters <- rev(unique(data_sf$quarter_year))
  
  # Initialize empty lists to store ppp and KDE objects
  region_quarters_conflict <- list()
  kde_plot_list <- list()
  
  # Step 2: Loop over each quarter and process conflict data
  for (quarter in quarters) {
    # Filter the data for the current quarter
    quarter_data <- data_sf %>%
      filter(quarter_year == quarter)
    
    # Extract coordinates for this quarter
    quarter_coords <- st_coordinates(quarter_data)
    
    # Create a ppp object for this quarter
    quarter_ppp <- ppp(
      x = quarter_coords[, 1],
      y = quarter_coords[, 2],
      window = region_window
    )
    
    # Step 3: Remove rejected points that fall outside the window
    valid_points <- inside.owin(quarter_ppp$x, quarter_ppp$y, region_window)
    quarter_ppp <- quarter_ppp[valid_points]
    
    # Step 4: Jitter duplicates to avoid overplotting
    is_duplicate <- duplicated(quarter_ppp)
    jittered_x <- quarter_ppp$x + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
    jittered_y <- quarter_ppp$y + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
    
    # Create a jittered ppp object with the new coordinates
    jittered_ppp <- ppp(
      x = jittered_x,
      y = jittered_y,
      window = region_window
    )
    
    # Step 5: Store the jittered ppp object for later use
    region_quarters_conflict[[quarter]] <- jittered_ppp
    
    # Step 6: Rescale to kilometers and calculate KDE with specified sigma type and kernel method
    jittered_ppp_km <- rescale(jittered_ppp, 1000, "km")
    
    # Use selected sigma type
    if (sigma_type == "scott") {
      sigma_value <- bw.scott(jittered_ppp_km)
      kde_quarter <- density(jittered_ppp_km, sigma = bw.scott(jittered_ppp_km), edge = TRUE, kernel = kernel_method)
    } else if (sigma_type == "diggle") {
      sigma_value <- bw.scott(jittered_ppp_km)
      kde_quarter <- density(jittered_ppp_km, sigma = bw.diggle(jittered_ppp_km), edge = TRUE, kernel = kernel_method)
    } else if (sigma_type == "ppl") {
      sigma_value <- bw.scott(jittered_ppp_km)
      kde_quarter <- density(jittered_ppp_km, sigma = bw.ppl(jittered_ppp_km), edge = TRUE, kernel = kernel_method)
    } else if (sigma_type == "cvl") {
      sigma_value <- bw.scott(jittered_ppp_km)
      kde_quarter <- density(jittered_ppp_km, sigma = bw.CvL(jittered_ppp_km), edge = TRUE, kernel = kernel_method)
    } else {
      stop("Invalid sigma_type specified.")
    }
    
    # Step 7: Store the KDE object for later use
    kde_plot_list[[quarter]] <- kde_quarter
    
    # Print a message after processing each quarter
    print(paste("KDE - Quarter:", quarter, "| Kernel:", kernel_method, "| Sigma:", sigma_type, " | Sigma Value: ", sigma_value))
  }
  
  # Return the processed ppp and KDE objects
  return(list(
    "ppp_list" = region_quarters_conflict,
    "kde_list" = kde_plot_list
  ))
}

7.2 Setting up a function to display the visual across quarters

Next, we will create a function to visualize the KDE layers for each quarter, enabling us to compare the intensity of conflict events over time. The function will:

  • Take the KDE layers from each quarter.

  • Display them side by side or as an animation to show how conflict density changes across quarters.

  • This visualization will help us observe trends, including whether conflicts are becoming more localized or dispersed over time.

# Visualization Function for Conflict Points and KDE
visualize_conflict_results <- function(results, region_name, sigma_type, kernel_method) {
  ppp_list <- results$ppp_list
  kde_list <- results$kde_list
  
  # Plotting Kernel Density Estimates (KDE) for each quarter
  kde_plots <- lapply(names(kde_list), function(quarter) {
    kde_data <- as.data.frame(as.im(kde_list[[quarter]]))
    
    ggplot() +
      geom_raster(data = kde_data, aes(x = x, y = y, fill = value), alpha = 0.8) +
      scale_fill_viridis_c(option = "inferno", name = "Density") +
      labs(title = paste("KDE for Civilian Conflict", region_name, "-", quarter)) +
      theme_void() +  # Removes axis, background, and grid
      theme(
        plot.title = element_text(hjust = 0.5, face = "bold", size = 10, margin = margin(t = 10, b = 10)),  # Centers and adds margin to the title
        plot.margin = margin(t = 5, r = 5, b = 5, l = 5)  # Adjusts the margin around the plot
      )
  })
  
  # Combine all plots into a grid layout with a main title
  kde_grid <- grid.arrange(grobs = kde_plots, ncol = 4, top = textGrob(
    paste(region_name, " Civilian Conflicts KDE By", sigma_type, "Using", kernel_method, "Kernel"),
    gp = gpar(fontface = "bold", fontsize = 16)
  ))
  
  return(kde_grid)
}

7.3 State’s Quarterly Kernel Density estimation.

For each state, we will calculate and visualize the Kernel Density Estimation (KDE) using two different kernels. We focus on (You can use any we are flexible in the function above!)

  • Kernel 1: Gaussian Kernel

  • Kernel 2: Epanechnikov Kernel

The Gaussian Kernel is smooth and gives a global view of the density, whereas the Epanechnikov Kernel is more efficient and provides sharper boundaries, making it easier to identify clusters or hotspots.

We will also compute the KDE with different bandwidths, which controls the smoothness of the density estimates. A larger bandwidth results in smoother KDE, while a smaller bandwidth captures more local variations.

7.3.1 Yangon

7.3.1.1 Kernel 1: Gaussian Kernel

Yangon_Results_Diggle_Gaussian <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "diggle",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
# Step 4: Visualize the results with a main title
visualizations <- visualize_conflict_results(Yangon_Results_Diggle_Gaussian, "Yangon", "Diggle", "Gaussian")

Yangon_Results_Scott_Gaussian <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "scott",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
# Step 4: Visualize the results with a main title
Yangon_Results_Scott_Gaussian_Visualisation <- visualize_conflict_results(Yangon_Results_Scott_Gaussian, "Yangon", "Scott", "Gaussian")

Yangon_Results_CvL_Gaussian <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "cvl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Yangon_Results_CvL_Gaussian_visualizations <- visualize_conflict_results(Yangon_Results_CvL_Gaussian, "Yangon", "CvL", "Gaussian")

Yangon_Results_PPL_Gaussian <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "ppl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
# Step 4: Visualize the results with a main title
Yangon_Results_PPL_Gaussian_visualizations <- visualize_conflict_results(Yangon_Results_PPL_Gaussian, "Yangon", "PPL", "Gaussian")

7.3.1.2 Kernel 2: Epanechnikov Kernel

Yangon_Results_Diggle_Epanechnikov <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "diggle",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Yangon_Results_Diggle_Epanechnikov_visualizations <- visualize_conflict_results(Yangon_Results_Diggle_Epanechnikov, "Yangon", "Diggle", "Epanechnikov")

Yangon_Results_Scott_Epanechnikov <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "scott",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Yangon_Results_Scot_Epanechnikov_visualizations <- visualize_conflict_results(Yangon_Results_Scott_Epanechnikov, "Yangon", "Scott", "Epanechnikov")

Yangon_Results_CvL_Epanechnikov <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "cvl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Yangon_Results_CvL_Epanechnikov_visualizations <- visualize_conflict_results(Yangon_Results_CvL_Epanechnikov, "Yangon", "CvL", "Epanechnikov")

Yangon_Results_PPL_Epanechnikov <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf,
  sigma_type = "ppl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Yangon_Results_PPL_Epanechnikov_visualizations <- visualize_conflict_results(Yangon_Results_PPL_Epanechnikov, "Yangon", "PPL", "Epanechnikov")

7.3.1.3 Analysis and Interpretation.

We used Scott’s bandwidth and the Epanechnikov kernel to model the spatial distribution of civilian conflicts in Yangon from 2021 to 2024. Scott’s method provides an optimal balance between bias and variance, ensuring smoother KDE without losing significant details. The Epanechnikov kernel minimizes error and focuses on areas with higher density, effectively capturing conflict hotspots while reducing the impact of outliers.

Key Observations:

  • Central Yangon Hotspot: Conflict density remains consistently high in central Yangon throughout the period.

  • Temporal Shifts: Conflict intensity fluctuates across quarters, peaking in certain periods (e.g., Q3-2022 and Q1-2024), with some spread into surrounding areas over time.

  • Density Spread: While central areas remain the most affected, the gradual spread suggests an intensification or expansion of conflicts.

In summary, Scott’s bandwidth and the Epanechnikov kernel provided a precise view of conflict clustering and temporal changes, revealing critical patterns in civilian-related conflict hotspots in Yangon.

7.3.2 Mandalay

7.3.2.1 Mandalay Kernel 1: Gaussian Kernel

Mandalay_Results_Diggle_Gaussian <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "diggle",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_Diggle_Gaussian_visualizations <- visualize_conflict_results(Mandalay_Results_Diggle_Gaussian, "Mandalay", "Diggle", "Gaussian")

Mandalay_Results_Scott_Gaussian <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "scott",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_Scott_Gaussian_visualizations <- visualize_conflict_results(Mandalay_Results_Scott_Gaussian, "Mandalay", "Scott", "Gaussian")

Mandalay_Results_CvL_Gaussian <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "cvl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_CvL_Gaussian_visualizations <- visualize_conflict_results(Mandalay_Results_CvL_Gaussian, "Mandalay", "CvL", "Gaussian")

Mandalay_Results_PPL_Gaussian <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "ppl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_PPL_Gaussian_visualizations <- visualize_conflict_results(Mandalay_Results_PPL_Gaussian, "Mandalay", "PPL", "Gaussian")

7.3.2.2 Mandalay Kernel 2: Epanechnikov Kernel

Mandalay_Results_Diggle_Epanechnikov <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "diggle",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_Diggle_Epanechnikov_visualizations <- visualize_conflict_results(Mandalay_Results_Diggle_Epanechnikov, "Mandalay ", "Diggle", "Epanechnikov")

Mandalay_Results_Scott_Epanechnikov <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "scott",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_Scott_Epanechnikov_visualizations <- visualize_conflict_results(Mandalay_Results_Scott_Epanechnikov, "Mandalay ", "Scott", "Epanechnikov")

Mandalay_Results_CvL_Epanechnikov <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay ",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "cvl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_CvL_Epanechnikov_visualizations <- visualize_conflict_results(Mandalay_Results_CvL_Epanechnikov, "Mandalay ", "CvL", "Epanechnikov")

Mandalay_Results_PPL_Epanechnikov <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf,
  sigma_type = "ppl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Mandalay_Results_PPL_Epanechnikov_visualizations <- visualize_conflict_results(Mandalay_Results_PPL_Epanechnikov, "Mandalay ", "PPL", "Epanechnikov")

7.3.2.3 Analysis and Interpretation

We applied Scott’s bandwidth and the Epanechnikov kernel for the Mandalay civilian conflict KDE, as it provides optimal smoothing and reduces variance. The kernel effectively focuses on dense areas, capturing conflict intensity over time without distortion from outliers.

Key Observations:

  • Central Mandalay Hotspot: Consistent conflict density is observed in central Mandalay across all quarters, with peaks in Q3-2022 and Q1-2023.

  • Spatiotemporal Patterns: Conflict density fluctuates quarter by quarter, with spread toward the north in several periods, indicating shifts in conflict zones.

  • Stability vs Spread: While central Mandalay remains a key hotspot, the spread of conflicts over time suggests areas of emerging concern.

This KDE approach reveals both stable hotspots and the evolving nature of conflicts in Mandalay, providing insights into the intensity and movement of conflict zones across different periods.

7.3.3 Sagaing

7.3.3.1 Sagaing Kernel 1: Gaussian Kernel

Sagaing_Results_Diggle_Gaussian <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "diggle",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_Diggle_Gaussian_visualizations <- visualize_conflict_results(Sagaing_Results_Diggle_Gaussian, "Sagaing", "Diggle", "Gaussian")

Sagaing_Results_Scott_Gaussian <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "scott",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_Scott_Gaussian_visualizations <- visualize_conflict_results(Sagaing_Results_Scott_Gaussian, "Sagaing", "Scott", "Gaussian")

Sagaing_Results_CVL_Gaussian <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "cvl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_CVL_Gaussian_visualizations <- visualize_conflict_results(Sagaing_Results_CVL_Gaussian, "Sagaing", "CVL", "Gaussian")

Sagaing_Results_PPL_Gaussian <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "ppl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "gaussian"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_PPL_Gaussian_visualizations <- visualize_conflict_results(Sagaing_Results_PPL_Gaussian, "Sagaing", "PPL", "Gaussian")

7.3.3.2 Sagaing Kernel 2: Epanechnikov Kernel

Sagaing_Results_Diggle_Epanechnikov <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "diggle",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_Diggle_Epanechnikov_visualizations <- visualize_conflict_results(Sagaing_Results_Diggle_Epanechnikov, "Sagaing", "Diggle", "Epanechnikov")

Sagaing_Results_Scott_Epanechnikov <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "scott",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_Scott_Epanechnikov_visualizations <- visualize_conflict_results(Sagaing_Results_Scott_Epanechnikov, "Sagaing", "Scott", "Epanechnikov")

Sagaing_Results_Cvl_Epanechnikov <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "cvl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_CvL_Epanechnikov_visualizations <- visualize_conflict_results(Sagaing_Results_Cvl_Epanechnikov, "Sagaing", "CvL", "Epanechnikov")

Sagaing_Results_PPL_Epanechnikov <- process_quarter_conflicts(
  region_sf = Sagaing_Sf,
  region_window = Sagaing_Owin,
  region_name = "Sagaing",
  data_sf = Sagaing_ACLED_Data_Sf,
  sigma_type = "ppl",  # Can change to "scott", "diggle", or "ppl"
  kernel_method = "epanechnikov"  # Can change to "epanechnikov", "quartic", or others
)
Sagaing_Results_PPL_Epanechnikov_visualizations <- visualize_conflict_results(Sagaing_Results_PPL_Epanechnikov, "Sagaing", "PPL", "Epanechnikov")

7.3.3.3 Analysis and Interpretation

In Sagaing, we again utilized Scott’s bandwidth and the Epanechnikov kernel for KDE to assess the spatiotemporal distribution of civilian conflict events. This method ensures that key dense areas are well-represented while minimizing the impact of noise from less significant areas.

Key Observations:

  • Consistent Hotspot: The southern region of Sagaing consistently shows high-density conflict areas throughout the quarters, especially peaking during Q3-2022 and Q1-2023.

  • Emerging Zones: Some quarters, such as Q1-2022 and Q3-2023, show increased conflict spread to the north, indicating a potential shift in conflict locations.

  • Stability and Peaks: The core conflict areas maintain their intensity, while certain periods experience an expansion of conflict areas, with a notable rise in density during the mid-2022 period.

These patterns highlight that conflict zones in Sagaing are relatively stable with periods of increased intensity, particularly in southern areas. This provides critical insights into when and where civilian conflict events have intensified.

8.0 Perform 2nd-Order Spatial Point Pattern Analysis

In this analysis, we aim to determine whether conflicts are clustered or randomly distributed within the study area. To achieve this, we use 2nd-order spatial point pattern analysis, which examines the relationships between individual conflict events at varying distances. This helps us detect patterns of clustering or dispersion across the state.

Why Use Only K and L Functions

We focus on the K Function and L Function because they provide the most straightforward and effective means of identifying clustering or dispersion. The K Function measures clustering at multiple scales, while the L Function simplifies interpretation by stabilizing variance, making it easier to identify spatial patterns in conflict events. These methods are ideal for capturing the spatial structure of conflicts across varying distances.

8.1 Yangon

For the Yangon conflict events, we’ll perform a detailed spatial point pattern analysis using the Clark-Evans test, Ripley’s K-function, and L-function. This will help us understand whether the conflict events are randomly distributed, clustered, or dispersed, and at what scales clustering occurs.

Null and Alternative Hypothesis:

  • Ho (Null Hypothesis): The distribution of civilian-related conflicts in Yangon is randomly distributed (CSR: Complete Spatial Randomness).

  • H1 (Alternative Hypothesis): The distribution of civilian-related conflicts in Yangon is not randomly distributed (it is clustered or dispersed).

8.1.1 Yangon’s Clark and Evan Test

# Perform the Clark-Evans test for clustering
Yangon_ClarksEvan <- clarkevans.test(Yangon_ACLED_ppp, 
                                      correction = "none", 
                                      clipregion = Yangon_Owin, 
                                      alternative = c("clustered"),
                                      nsim = 30)

# Print Clark-Evans p-value and R-statistic
print(Yangon_ClarksEvan)

    Clark-Evans test
    No edge correction
    Z-test

data:  Yangon_ACLED_ppp
R = 0.11849, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

8.1.2 Yangon’s K-Function Method

K_Yangon <- Kest(Yangon_ACLED_ppp, correction = "Ripley")
# Step 2: Plot the K-function, showing K(d) - r
plot(K_Yangon, . -r ~ r, ylab = "K(d) - r", xlab = "d(KM)", main = "Ripley's K-Function for Yangon Civilian Related Conflicts")

K_Yangon_CSR <- envelope(Yangon_ACLED_ppp, Kest, nsim = 30, rank = 1, global = TRUE)
plot(K_Yangon_CSR, . - r ~ r, xlab = "d(KM", ylab = "K(d) - r", main = "Envelope for K-Function (CSR) - Yangon Civilian Related Conflicts")

8.1.2.2 Yangon’s K-Function Method Observation Table

# Step 1: Define key distances (e.g., 20000m, 30000m, 40000m)
key_distances <- c(5, 10, 15,20,25)

# Step 2: Extract observed K-function values and CSR bounds at distances closest to key distances
closest_indices <- sapply(key_distances, function(d) which.min(abs(K_Yangon_CSR$r - d)))

# Step 3: Create a table summarizing observed and CSR envelope bounds at the closest distances
observed_table <- data.frame(
  Distance = K_Yangon_CSR$r[closest_indices],       # Actual distances used in the analysis
  K_Obs = K_Yangon_CSR$obs[closest_indices],        # Observed K-function values
  K_Lo = K_Yangon_CSR$lo[closest_indices],          # Lower bound of CSR envelope
  K_Hi = K_Yangon_CSR$hi[closest_indices]           # Upper bound of CSR envelope
)

# Display the observed table
print(observed_table)
  Distance     K_Obs       K_Lo      K_Hi
1  4.99309  1822.917   35.71503  120.9307
2  9.98618  5226.240  270.68366  355.8994
3 14.97927  8118.862  662.29805  747.5137
4 20.02663  9683.831 1217.37812 1302.5938
5 25.01972 10687.263 1923.98668 2009.2024

8.1.2 Yangon’s L Functions

L_Yangon <- Lest(Yangon_ACLED_ppp, correction = "Ripley")
plot(L_Yangon, . -r ~ r, ylab = "L(d) - r", xlab = "d(KM)", main = "Ripley's L-Function for Yangon Civilian Related Conflicts")

L_Yangon_CSR <- envelope(Yangon_ACLED_ppp, Lest, nsim = 30, rank = 1, global = TRUE, savefuns = TRUE)
plot(L_Yangon_CSR, . - r ~ r, xlab = "d(km)", ylab = "L(d) - r", main = "Envelope for L-Function (CSR) - Yangon Civilian Related Conflicts")

8.1.2.2 Yangon’s L Functions Observation Table

key_distances <- c(5, 10, 15,20,25)
# Step 2: Extract observed K-function values and CSR bounds at distances closest to key distances
closest_indices <- sapply(key_distances, function(d) which.min(abs(L_Yangon_CSR$r - d)))
# Step 3: Create a table summarizing observed and CSR envelope bounds at the closest distances
observed_table <- data.frame(
  Distance = L_Yangon_CSR$r[closest_indices],       # Actual distances used in the analysis
  L_Obs = L_Yangon_CSR$obs[closest_indices],        # Observed K-function values
  L_Lo = L_Yangon_CSR$lo[closest_indices],          # Lower bound of CSR envelope
  L_Hi = L_Yangon_CSR$hi[closest_indices]           # Upper bound of CSR envelope
)
# Display the observed table
print(observed_table)
  Distance    L_Obs      L_Lo      L_Hi
1  4.99309 24.08843  4.836161  5.150019
2  9.98618 40.78681  9.829251 10.143108
3 14.97927 50.83615 14.822341 15.136198
4 20.02663 55.51990 19.869703 20.183561
5 25.01972 58.32548 24.862793 25.176651

8.1.4 Yangons- Analysis and Interpretation

Based on the analysis of the Clark-Evans test, K-function, and L-function, we can draw the following conclusions:

  1. Clark-Evans Test:

    • R = 0.099924, p-value < 2.2e-16 . The R-value is less than 1, and the p-value is very small, which means the civilian related conflict events in Yangon are not randomly distributed and are significantly clustered.

  2. K-function Analysis:

    • The observed K-function consistently exceeds the upper bound of the CSR envelope across all distances, indicating that the clustering is prominent over a wide range of distances (up to 25km). The clustering becomes more pronounced as the distance increases.

  3. L-function Analysis:

    x

    • The observed L-function also exceeds the CSR envelope’s upper bounds, confirming the presence of clustering at multiple spatial scales. This supports the findings from the K-function.

    1. Observed Tables:

      Distance K_Obs K_Lo K_Hi L_OBV L_Lo L_Hi
      4.99309 1822.917 35.71503 120.9307 24.0883 4.836161 5.150019
      9.98618 5226.240 270.68366 355.8994 40.78681 9.829251 10.14311
      14.97927 8118.682 662.29806 747.5137 50.83615 14.82234 15.1362
      20.02663 9683.831 1217.37812 1302.5938 55.5199 19.8697 20.18356
      25.01972 10687.263 1923.98668 2009.202 53.32548 24.86279 25.17665
      • The observed K-function and L-function tables confirm that the observed values are significantly higher than the expected values under CSR, particularly at distances of 5km, 10km, 15km, and beyond, suggesting clustering at larger spatial scales.

We reject the null hypothesis of Complete Spatial Randomness (CSR) for Yangon civilian-related conflicts. The Clark-Evans test shows significant clustering, and both the K-function and L-function support this finding by indicating clustering over a broad range of distances. The observed tables further confirm this, as the observed values consistently exceed the CSR bounds. Therefore, the civilian conflict events in Yangon are significantly clustered rather than randomly distributed.

8.2 Mandalay

For the Mandalay conflict events, we’ll perform a detailed spatial point pattern analysis using the Clark-Evans test, Ripley’s K-function, and L-function. This will help us understand whether the conflict events are randomly distributed, clustered, or dispersed, and at what scales clustering occurs.

Null and Alternative Hypothesis:

  • Ho (Null Hypothesis): The distribution of civilian-related conflicts in Mandalay is randomly distributed (CSR: Complete Spatial Randomness).

  • H1 (Alternative Hypothesis): The distribution of civilian-related conflicts in Mandalay is not randomly distributed (it is clustered or dispersed).

8.2.1 Clarks Evan Test

# Perform the Clark-Evans test for clustering
Mandalay_ClarksEvan <- clarkevans.test(Mandalay_ACLED_ppp, 
                              correction = "none", 
                              clipregion = Mandalay_Owin, 
                              alternative = c("clustered"),
                              nsim = 30)
# Print Clark-Evans p-value and R-statistic
print(Mandalay_ClarksEvan)

    Clark-Evans test
    No edge correction
    Z-test

data:  Mandalay_ACLED_ppp
R = 0.2192, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

8.2.2 K Function

K_Mandalay <- Kest(Mandalay_ACLED_ppp, correction = "Ripley")
plot(K_Mandalay, . -r ~ r, ylab = "K(d) - r", xlab = "d(KM)", main = "Ripley's K-Function for Mandalay Civilian Related Conflicts")

K_Mandalay_CSR <- envelope(Mandalay_ACLED_ppp, Kest, nsim = 30, rank = 1, global = TRUE)
plot(K_Mandalay_CSR, . -r ~ r, ylab = "K(d) - r", xlab = "d(Km)", main = "Ripley's K-Function for Mandalay Civilian Related Conflicts")

8.2.2.2 Observed Table (K-Function at Key Distances)

# Step 1: Define key distances (e.g., 20000m, 30000m, 40000m)
key_distances <- c(10,20,30,40)

# Step 2: Extract observed K-function values and CSR bounds at distances closest to key distances
closest_indices <- sapply(key_distances, function(d) which.min(abs(K_Mandalay_CSR$r - d)))

# Step 3: Create a table summarizing observed and CSR envelope bounds at the closest distances
observed_table <- data.frame(
  Distance = K_Mandalay_CSR$r[closest_indices],       # Actual distances used in the analysis
  K_Obs = K_Mandalay_CSR$obs[closest_indices],        # Observed K-function values
  K_Lo = K_Mandalay_CSR$lo[closest_indices],          # Lower bound of CSR envelope
  K_Hi = K_Mandalay_CSR$hi[closest_indices]           # Upper bound of CSR envelope
)

# Display the observed table
print(observed_table)
   Distance     K_Obs      K_Lo      K_Hi
1  9.962643  7305.334  260.9636  362.6693
2 20.016687 10852.795 1207.8820 1309.5878
3 29.979331 15847.857 2772.6858 2874.3915
4 40.033374 21114.336 4984.0867 5085.7925

8.2.3 L Function

L_Mandalay <- Lest(Mandalay_ACLED_ppp, correction = "Ripley")
plot(L_Mandalay, . -r ~ r, ylab = "L(d) - r", xlab = "d(KM)", main = "Ripley's L-Function for Mandalay Civilian Related Conflicts")

L_Mandalay_CSR <- envelope(Mandalay_ACLED_ppp, Lest, nsim = 30, rank = 1, global = TRUE, savefuns = TRUE)
plot(L_Mandalay_CSR, . - r ~ r, xlab = "d(KM)", ylab = "L(d) - r", main = "Envelope for L-Function (CSR) - Mandalay Civilian Related Conflicts")

8.2.3.2 Observed Table (L-Function at Key Distances)

key_distances <- c(10,20,30,40)

# Step 2: Extract observed K-function values and CSR bounds at distances closest to key distances
closest_indices <- sapply(key_distances, function(d) which.min(abs(L_Mandalay_CSR$r - d)))

# Step 3: Create a table summarizing observed and CSR envelope bounds at the closest distances
observed_table <- data.frame(
  Distance = L_Mandalay_CSR$r[closest_indices],       # Actual distances used in the analysis
  L_Obs = L_Mandalay_CSR$obs[closest_indices],        # Observed K-function values
  L_Lo = L_Mandalay_CSR$lo[closest_indices],          # Lower bound of CSR envelope
  L_Hi = L_Mandalay_CSR$hi[closest_indices]           # Upper bound of CSR envelope
)

# Display the observed table
print(observed_table)
   Distance    L_Obs      L_Lo     L_Hi
1  9.962643 48.22199  9.689332 10.23595
2 20.016687 58.77544 19.743376 20.29000
3 29.979331 71.02485 29.706020 30.25264
4 40.033374 81.98111 39.760063 40.30669

8.2.4 Mandalay- Analysis and Interpretation

Based on the analysis of the Clark-Evans test, K-function, and L-function, we can draw the following conclusions:

  1. Clark-Evans Test:

    • R = 0.14898, p-value < 2.2e-16. The R-value is less than 1, and the p-value is very small, which means the civilian related conflict events in Mandalay are not randomly distributed and are significantly clustered.

  2. K-function Analysis:

    • The observed K-function consistently exceeds the upper bound of the CSR envelope across all distances, indicating that the clustering is prominent over a wide range of distances (up to 70km). The clustering becomes more pronounced as the distance increases.

  3. L-function Analysis:

    • The observed L-function also exceeds the CSR envelope’s upper bounds, confirming the presence of clustering at multiple spatial scales. This supports the findings from the K-function.

  4. Observed Tables:

Distance K_Obs K_Lo K_Hi L_Obs L_Lo L_Hi
9.962643 7305.334 260.9636 362.6693 48.22199 9.689332 10.23595
20.016687 10852.795 1207.8820 1309.5878 58.77544 19.743376 20.2900
29.979331 15847.857 2772.6858 2874.3915 71.02485 29.706020 30.25264
40.033374 21114.336 4984.0867 5085.7925 81.9811 39.760063 40.30
  • The observed K-function and L-function tables confirm that the observed values are significantly higher than the expected values under CSR, particularly at distances of 20km, 30km, 40km, and beyond, suggesting clustering at larger spatial scales.

We reject the null hypothesis of Complete Spatial Randomness (CSR) for Mandalay civilian-related conflicts. The Clark-Evans test shows significant clustering, and both the K-function and L-function support this finding by indicating clustering over a broad range of distances. The observed tables further confirm this, as the observed values consistently exceed the CSR bounds. Therefore, the civilian conflict events in Mandalay are significantly clustered rather than randomly distributed.

8.3 Sagaing

For the Sagaing conflict events, we’ll perform a detailed spatial point pattern analysis using the Clark-Evans test, Ripley’s K-function, and L-function. This will help us understand whether the conflict events are randomly distributed, clustered, or dispersed, and at what scales clustering occurs.

Null and Alternative Hypothesis:

  • Ho (Null Hypothesis): The distribution of civilian-related conflicts in Sagaing is randomly distributed (CSR: Complete Spatial Randomness).

  • H1 (Alternative Hypothesis): The distribution of civilian-related conflicts in Sagaing is not randomly distributed (it is clustered or dispersed).

8.3.1 Clarks Evan Test

# Perform the Clark-Evans test for clustering
Sagaing_ClarksEvan <- clarkevans.test(Sagaing_ACLED_ppp, 
                              correction = "none", 
                              clipregion = Sagaing_Owin, 
                              alternative = c("clustered"),
                              nsim = 30)
# Print Clark-Evans p-value and R-statistic
print(Sagaing_ClarksEvan)

    Clark-Evans test
    No edge correction
    Z-test

data:  Sagaing_ACLED_ppp
R = 0.1827, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
K_Sagaing <- Kest(Sagaing_ACLED_ppp, correction = "Ripley")
plot(K_Sagaing, . -r ~ r, ylab = "K(d) - r", xlab = "d(KM)", main = "Ripley's K-Function for Sagaing Civilian Related Conflicts")

K_Sagaing_CSR <- envelope(Sagaing_ACLED_ppp, Kest, nsim = 30, rank = 1, global = TRUE)
plot(K_Sagaing_CSR, . - r ~ r, xlab = "d(KM)", ylab = "K(d) - r", main = "Envelope for K-Function (CSR) - Sagaing Conflicts")

8.3.2.2 Observed Table (K-Function at Key Distances)

# Step 1: Define key distances (e.g., 20000m, 30000m, 40000m)
key_distances <- c(20,30,40,50,60,70)

# Step 2: Extract observed K-function values and CSR bounds at distances closest to key distances
closest_indices <- sapply(key_distances, function(d) which.min(abs(K_Sagaing_CSR$r - d)))

# Step 3: Create a table summarizing observed and CSR envelope bounds at the closest distances
observed_table <- data.frame(
  Distance = K_Sagaing_CSR$r[closest_indices],       # Actual distances used in the analysis
  K_Obs = K_Sagaing_CSR$obs[closest_indices],        # Observed K-function values
  K_Lo = K_Sagaing_CSR$lo[closest_indices],          # Lower bound of CSR envelope
  K_Hi = K_Sagaing_CSR$hi[closest_indices]           # Upper bound of CSR envelope
)

# Display the observed table
print(observed_table)
  Distance     K_Obs      K_Lo      K_Hi
1 19.99394  7021.421  1186.336  1325.415
2 29.99091 12903.634  2756.181  2895.260
3 39.98788 19530.590  4953.964  5093.043
4 49.98485 26873.240  7779.685  7918.764
5 59.98182 34286.168 11233.343 11372.422
6 68.24599 40939.246 14562.475 14701.554

8.3.3 L-Function Method

# Step 1: Calculate the L-function (Ripley's L) for Sagaing conflicts
L_Sagaing <- Lest(Sagaing_ACLED_ppp, correction = "Ripley")
# Step 2: Plot the L-function, showing L(d) - r
plot(L_Sagaing, . - r ~ r, ylab = "L(d) - r", xlab = "d(km)", main = "Ripley's L-Function for Sagaing Civilian Related Conflicts")

L_Sagaing_CSR <- envelope(Sagaing_ACLED_ppp, Lest, nsim = 30, rank = 1, global = TRUE)
plot(L_Sagaing_CSR, . - r ~ r, xlab = "d(KM)", ylab = "K(d) - r", main = "Envelope for L-Function (CSR) - Sagaing Civilian Related Conflicts")

8.3.3.3 Observed Table (L-Function at Key Distances)

key_distances <- c(20,30,40,50,60,70)

# Step 2: Extract observed K-function values and CSR bounds at distances closest to key distances
closest_indices <- sapply(key_distances, function(d) which.min(abs(L_Sagaing_CSR$r - d)))

# Step 3: Create a table summarizing observed and CSR envelope bounds at the closest distances
observed_table <- data.frame(
  Distance = L_Sagaing_CSR$r[closest_indices],       # Actual distances used in the analysis
  L_Obs = L_Sagaing_CSR$obs[closest_indices],        # Observed K-function values
  L_Lo = L_Sagaing_CSR$lo[closest_indices],          # Lower bound of CSR envelope
  L_Hi = L_Sagaing_CSR$hi[closest_indices]           # Upper bound of CSR envelope
)

# Display the observed table
print(observed_table)
  Distance     L_Obs     L_Lo     L_Hi
1 19.99394  47.27566 19.81300 20.17488
2 29.99091  64.08864 29.80997 30.17185
3 39.98788  78.84656 39.80694 40.16882
4 49.98485  92.48793 49.80391 50.16579
5 59.98182 104.46830 59.80088 60.16277
6 68.24599 114.15501 68.06505 68.42693

8.3.4 Sagaing - Analysis and Interpretation

Based on the analysis of the Clark-Evans test, K-function, and L-function, we can draw the following conclusions:

  1. Clark-Evans Test:

    • R = 0.1485, p-value < 2.2e-16 . The R-value is less than 1, and the p-value is very small, which means the conflict events in Sagaing are not randomly distributed and are significantly clustered.

  2. K-function Analysis:

    • The observed K-function consistently exceeds the upper bound of the CSR envelope across all distances, indicating that the clustering is prominent over a wide range of distances (up to 70km). The clustering becomes more pronounced as the distance increases.

  3. L-function Analysis:

    • The observed L-function also exceeds the CSR envelope’s upper bounds, confirming the presence of clustering at multiple spatial scales. This supports the findings from the K-function.

  4. Observed Tables:

Distance K_Obs K_Lo K_Hi L_Obs L_Lo L_Hi
19.99394 7021.421 1186.336 1325.415 19.81300 20.17488 20.17488
29.99091 12903.634 2756.181 2895.260 64.08864 29.80997 30.17185
39.98788 19530.590 4963.964 5093.043 78.84656 39.80694 40.16882
49.98485 26873.240 7797.685 7916.764 92.48793 49.80391 50.16579
59.98182 34286.168 11233.343 11372.422 104.46830 59.80088 60.16277
68.24599 40939.246 1456,475 14701.554 114.15501 68.06505 68.4269
  • The observed K-function and L-function tables confirm that the observed values are significantly higher than the expected values under CSR, particularly at distances of 20km, 30km, 40km, and beyond, suggesting clustering at larger spatial scales.

We reject the null hypothesis of Complete Spatial Randomness (CSR) for Sagaing civilian-related conflicts. The Clark-Evans test shows significant clustering, and both the K-function and L-function support this finding by indicating clustering over a broad range of distances. The observed tables further confirm this, as the observed values consistently exceed the CSR bounds. Therefore, the conflict events in Sagaing are significantly clustered rather than randomly distributed.

9.0 Derive Quarterly Spatio-Temporal KDE Layers

In this section, we will generate Quarterly Spatio-Temporal Kernel Density Estimation (KDE) layers to analyze how conflict intensity evolves across Yangon, Mandalay, and Sagaing. By applying KDE across different quarters, we aim to capture the spatial distribution and temporal shifts in conflict hotspots. This will help us track changes in the concentration and spread of events, providing insights into regional conflict dynamics over time.

9.1 Steps for code Breakdown.

  • Step 1: Convert data to sf object to ensure compatibility for spatial analysis.

  • Step 2: Extract longitude and latitude from the geometry column for plotting and analysis.

  • Step 3: Create numeric and factor representations of quarters for time-based KDE and distinct plotting.

  • Step 4: Generate a point pattern object (ppp) for spatial analysis with quarter factors as marks.

  • Step 5: Assign circle sizes to each quarter for visual differentiation.

  • Step 6: Plot the base map (Yangon) without data points.

  • Step 7: Add points with sizes based on quarters to the map for visualization.

  • Step 8: Add a legend to explain the circle sizes and quarters.

  • Step 9: Create a numeric point pattern for KDE based on quarter information.

  • Step 10: Perform KDE analysis on the spatio-temporal point pattern.

  • Step 11: Define the quarters to be plotted for KDE analysis.

  • Step 12: Set up a 4x4 plotting grid to display multiple KDE plots.

  • Step 13: Loop through quarters and plot KDE layers for each, visualizing conflict intensity over time.

9.2 Yangon Quarter Year Spatio-Temporal KDE

# Step 1: Ensure Yangon_QuarterYear_Sf is an sf object (if not already)
Yangon_QuarterYear_Sf <- st_as_sf(Yangon_ACLED_Data_Sf)

# Step 2: Extract the coordinates (longitude and latitude) from the geometry column
coords <- st_coordinates(Yangon_QuarterYear_Sf)

# Step 3: Convert quarter_numeric to factor for distinct quarters (for plotting the circles)
Yangon_QuarterYear_Sf <- Yangon_QuarterYear_Sf %>%
  mutate(
    year = as.numeric(sub(".*-", "", quarter_year)),                  # Extract the year (e.g., 2023)
    quarter = as.numeric(sub("Q", "", sub("-.*", "", quarter_year))), # Extract the quarter (e.g., Q1 -> 1)
    
    # Continuous representation of time (factor for distinct quarters)
    quarter_numeric_factor = as.factor(year * 10 + quarter),  # Convert to factor for plotting the circles
    quarter_numeric = year * 10 + quarter  # Keep as numeric for time-based KDE analysis
  )

# Step 4: Create the point pattern object using ppp() with factor marks (for distinct quarter plotting)
Yangon_QuarterYear_PPP_Factor <- ppp(
  x = coords[, 1],  # Longitude (x-coordinates)
  y = coords[, 2],  # Latitude (y-coordinates)
  window = Yangon_Owin,  # Spatial window (Yangon_Owin)
  marks = Yangon_QuarterYear_Sf$quarter_numeric_factor  # Use factor for visualizing distinct quarters
)


# Step 5: Assign 14 unique sizes to the 14 unique quarter_numeric values
# Create 14 different circle sizes, e.g., from 1 to 3 in size
unique_quarters <- levels(Yangon_QuarterYear_Sf$quarter_numeric_factor)  # Get the unique quarter levels (14 levels)
circle_sizes <- seq(1, 3, length.out = length(unique_quarters))  # Generate 14 sizes from 1 to 3

# Map each unique quarter_numeric level to a circle size
circle_size_map <- setNames(circle_sizes, unique_quarters)  # Create a named vector to map quarter to size
# Apply the circle size mapping to each point based on its quarter_numeric_factor
point_circle_sizes <- circle_size_map[Yangon_QuarterYear_Sf$quarter_numeric_factor]

# Step 6: Plot the base map (Yangon_Owin) without points
plot(Yangon_Owin, main = "Yangon Quarter-Year Owin with Unique Circle Sizes")

# Step 7: Add the points with the unique circle sizes based on quarter_numeric
points(coords[, 1], coords[, 2], cex = point_circle_sizes, pch = 16, col = "black")  # Add circles on top of base map

# Step 8: Add a custom legend to show the quarter_numeric values and corresponding circle sizes
legend("topright", legend = unique_quarters, 
       pch = 16, 
       pt.cex = circle_sizes,  # Show unique circle sizes in the legend
       col = "black", 
       title = "Quarters (Circle Size)")

Yangon_QuarterYear_PPP_Numeric <- ppp(
  x = coords[, 1],  # Longitude (x-coordinates)
  y = coords[, 2],  # Latitude (y-coordinates)
  window = Yangon_Owin,  # Spatial window (Yangon_Owin)
  marks = Yangon_QuarterYear_Sf$quarter_numeric  # Use numeric for KDE analysis (time-based)
)
Yangon_QuarterYear_PPP_Numeric <- rescale(Yangon_QuarterYear_PPP_Numeric, 1000, "km")


Yangon_QuarterYear_KDE <- spattemp.density(Yangon_QuarterYear_PPP_Numeric)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
multi_color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))


# Define the quarters for which you want to plot the KDE
tims <- c(20211, 20212, 20213, 20214, 20221, 20222, 20223, 20224, 
          20231, 20232, 20233, 20234, 20241, 20242)

#Set up the plotting window to display multiple plots in a grid (4 rows, 4 columns)
par(mfrow=c(4, 4), mar=c(2, 2, 2, 2), oma=c(0, 0, 2, 0))  # Adjust margins and add outer margin for title

# Loop through the 'tims' vector and plot the KDE for each time point (from left to right)
for(i in tims) { 
    plot(Yangon_QuarterYear_KDE, i, 
         override.par=FALSE,  # Keep the graphical parameters unchanged
         fix.range=TRUE,      # Fix the range of the KDE
         main=paste("Quarter", i),  # Title for each plot
         ribside = "right", 
         col = multi_color_palette(100))  # Use 'inferno' color gradient
         
    # Remove the x and y axis labels for a cleaner look
    axis(1, labels=FALSE)
    axis(2, labels=FALSE)
}

#dd a unified title across all plots (e.g., the duration of analysis)
mtext("Quarterly Spatio-Temporal KDE for Civilian Conflict Intensity in Yangon", 
      outer=TRUE, cex=1.5, font=2)

multi_color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))

# Save the animation as a GIF
saveGIF({
  # Loop over the valid quarter times
  for(i in tims){ 
    suppressWarnings({
      plot(Yangon_QuarterYear_KDE, i, 
           override.par=FALSE,  # Keep graphical parameters unchanged
           fix.range=TRUE,      # Fix the range of the KDE
           main=paste("Mandalay Civilian Related Conflict KDE at Quarter", i),  # Title for each plot
           ribside = "right",  # Legend on the right
           col = multi_color_palette(100))  # Apply multi-color gradient
    })
  }
}, movie.name = "yangonkde_animation.gif", interval = 0.5, ani.width = 800, ani.height = 800)

In Yangon, the spatio-temporal KDE analysis shows a reduction in overall conflict intensity from 2021 to 2024. However, it is important to note that while the intensity has decreased, clustering is still evident within the same few areas.

Key Observations:

  • 2021 to 2022: Widespread high-intensity conflict areas in central Yangon are prominent in 2021 and early 2022, indicated by the larger red and orange zones.

  • 2023 to 2024: Although the intensity diminishes over time, the clustering remains persistent in specific areas. The high-intensity zones shrink, but they consistently appear in similar locations, suggesting that while the scale of conflict has reduced, these areas remain focal points for civilian conflict.

The persistent clustering, even as intensity diminishes, indicates that certain regions continue to experience conflicts, requiring targeted interventions in those zones.

9.3 Mandalay Quarter Year Spatio-Temporal KDE

# Step 1: Ensure Mandalay_QuarterYear_Sf is an sf object (if not already)
Mandalay_QuarterYear_Sf <- st_as_sf(Mandalay_ACLED_Data_Sf)

# Step 2: Extract the coordinates (longitude and latitude) from the geometry column
coords <- st_coordinates(Mandalay_QuarterYear_Sf)

# Step 3: Convert quarter_numeric to factor for distinct quarters (for plotting the circles)
Mandalay_QuarterYear_Sf <- Mandalay_QuarterYear_Sf %>%
  mutate(
    year = as.numeric(sub(".*-", "", quarter_year)),                  # Extract the year (e.g., 2023)
    quarter = as.numeric(sub("Q", "", sub("-.*", "", quarter_year))), # Extract the quarter (e.g., Q1 -> 1)
    
    # Continuous representation of time (factor for distinct quarters)
    quarter_numeric_factor = as.factor(year * 10 + quarter),  # Convert to factor for plotting the circles
    quarter_numeric = year * 10 + quarter  # Keep as numeric for time-based KDE analysis
  )

# Step 4: Create the point pattern object using ppp() with factor marks (for distinct quarter plotting)
Mandalay_QuarterYear_PPP_Factor <- ppp(
  x = coords[, 1],  # Longitude (x-coordinates)
  y = coords[, 2],  # Latitude (y-coordinates)
  window = Mandalay_Owin,  # Spatial window (Mandalay_Owin)
  marks = Mandalay_QuarterYear_Sf$quarter_numeric_factor  # Use factor for visualizing distinct quarters
)


# Step 5: Assign 14 unique sizes to the 14 unique quarter_numeric values
# Create 14 different circle sizes, e.g., from 1 to 3 in size
unique_quarters <- levels(Mandalay_QuarterYear_Sf$quarter_numeric_factor)  # Get the unique quarter levels (14 levels)
circle_sizes <- seq(1, 3, length.out = length(unique_quarters))  # Generate 14 sizes from 1 to 3

# Map each unique quarter_numeric level to a circle size
circle_size_map <- setNames(circle_sizes, unique_quarters)  # Create a named vector to map quarter to size
# Apply the circle size mapping to each point based on its quarter_numeric_factor
point_circle_sizes <- circle_size_map[Mandalay_QuarterYear_Sf$quarter_numeric_factor]

# Step 6: Plot the base map (Mandalay_Owin) without points
plot(Mandalay_Owin, main = "Mandalay Quarter-Year Owin with Unique Circle Sizes")

# Step 7: Add the points with the unique circle sizes based on quarter_numeric
points(coords[, 1], coords[, 2], cex = point_circle_sizes, pch = 16, col = "black")  # Add circles on top of base map

# Step 8: Add a custom legend to show the quarter_numeric values and corresponding circle sizes
legend("topright", legend = unique_quarters, 
       pch = 16, 
       pt.cex = circle_sizes,  # Show unique circle sizes in the legend
       col = "black", 
       title = "Quarters (Circle Size)")

Mandalay_QuarterYear_PPP_Numeric <- ppp(
  x = coords[, 1],  # Longitude (x-coordinates)
  y = coords[, 2],  # Latitude (y-coordinates)
  window = Mandalay_Owin,  # Spatial window (Mandalay_Owin)
  marks = Mandalay_QuarterYear_Sf$quarter_numeric  # Use numeric for KDE analysis (time-based)
)

Mandalay_QuarterYear_KDE <- spattemp.density(Mandalay_QuarterYear_PPP_Numeric)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
multi_color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))


# Define the quarters for which you want to plot the KDE
tims <- c(20211, 20212, 20213, 20214, 20221, 20222, 20223, 20224, 
          20231, 20232, 20233, 20234, 20241, 20242)

#Set up the plotting window to display multiple plots in a grid (4 rows, 4 columns)
par(mfrow=c(4, 4), mar=c(2, 2, 2, 2), oma=c(0, 0, 2, 0))  # Adjust margins and add outer margin for title

# Loop through the 'tims' vector and plot the KDE for each time point (from left to right)
for(i in tims) { 
  plot(Mandalay_QuarterYear_KDE, i, 
       override.par=FALSE,  # Keep the graphical parameters unchanged
       fix.range=TRUE,      # Fix the range of the KDE
       main=paste("Quarter", i),  # Title for each plot
       ribside = "right", 
       col = multi_color_palette(100))  # Use 'inferno' color gradient
  
  # Remove the x and y axis labels for a cleaner look
  axis(1, labels=FALSE)
  axis(2, labels=FALSE)
}

#dd a unified title across all plots (e.g., the duration of analysis)
mtext("Quarterly Spatio-Temporal KDE for Civilian Conflict Intensity in Mandalay", 
      outer=TRUE, cex=1.5, font=2)

multi_color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))

# Save the animation as a GIF
saveGIF({
  # Loop over the valid quarter times
  for(i in tims){ 
    suppressWarnings({
      plot(Mandalay_QuarterYear_KDE, i, 
           override.par=FALSE,  # Keep graphical parameters unchanged
           fix.range=TRUE,      # Fix the range of the KDE
           main=paste("Mandalay Civilian Related Conflict KDE at Quarter", i),  # Title for each plot
           ribside = "right",  # Legend on the right
           col = multi_color_palette(100))  # Apply multi-color gradient
    })
  }
}, movie.name = "Mandalaykde_animation.gif", interval = 0.5, ani.width = 800, ani.height = 800)

In Mandalay, the spatio-temporal KDE analysis similarly shows a general reduction in conflict intensity over time, from 2021 to 2024. However, like in Yangon, clustering remains within the same key areas across the quarters.

Key Observations:

  • 2021 to 2022: High-intensity conflicts are concentrated in central and northern Mandalay, as shown by the red and orange clusters, particularly in early 2021 and 2022.

  • 2023 to 2024: The intensity of the conflicts decreases over time, but these clusters persist in similar regions, particularly in central Mandalay. This indicates that although conflict intensity is declining, these areas continue to be hotspots for civilian conflict.

Thus, the persistence of clustering in the same areas despite decreasing conflict intensity highlights the need for focused intervention strategies in these conflict-prone regions.

9.4 Sagaing Quarter Year Spatio-Temporal KDE

# Step 1: Ensure Sagaing_QuarterYear_Sf is an sf object (if not already)
Sagaing_QuarterYear_Sf <- st_as_sf(Sagaing_ACLED_Data_Sf)

# Step 2: Extract the coordinates (longitude and latitude) from the geometry column
coords <- st_coordinates(Sagaing_QuarterYear_Sf)

# Step 3: Convert quarter_numeric to factor for distinct quarters (for plotting the circles)
Sagaing_QuarterYear_Sf <- Sagaing_QuarterYear_Sf %>%
  mutate(
    year = as.numeric(sub(".*-", "", quarter_year)),                  # Extract the year (e.g., 2023)
    quarter = as.numeric(sub("Q", "", sub("-.*", "", quarter_year))), # Extract the quarter (e.g., Q1 -> 1)
    
    # Continuous representation of time (factor for distinct quarters)
    quarter_numeric_factor = as.factor(year * 10 + quarter),  # Convert to factor for plotting the circles
    quarter_numeric = year * 10 + quarter  # Keep as numeric for time-based KDE analysis
  )

# Step 4: Create the point pattern object using ppp() with factor marks (for distinct quarter plotting)
Sagaing_QuarterYear_PPP_Factor <- ppp(
  x = coords[, 1],  # Longitude (x-coordinates)
  y = coords[, 2],  # Latitude (y-coordinates)
  window = Sagaing_Owin,  # Spatial window (Sagaing_Owin)
  marks = Sagaing_QuarterYear_Sf$quarter_numeric_factor  # Use factor for visualizing distinct quarters
)


# Step 5: Assign 14 unique sizes to the 14 unique quarter_numeric values
# Create 14 different circle sizes, e.g., from 1 to 3 in size
unique_quarters <- levels(Sagaing_QuarterYear_Sf$quarter_numeric_factor)  # Get the unique quarter levels (14 levels)
circle_sizes <- seq(1, 3, length.out = length(unique_quarters))  # Generate 14 sizes from 1 to 3

# Map each unique quarter_numeric level to a circle size
circle_size_map <- setNames(circle_sizes, unique_quarters)  # Create a named vector to map quarter to size
# Apply the circle size mapping to each point based on its quarter_numeric_factor
point_circle_sizes <- circle_size_map[Sagaing_QuarterYear_Sf$quarter_numeric_factor]

# Step 6: Plot the base map (Sagaing_Owin) without points
plot(Sagaing_Owin, main = "Sagaing Quarter-Year Owin with Unique Circle Sizes")

# Step 7: Add the points with the unique circle sizes based on quarter_numeric
points(coords[, 1], coords[, 2], cex = point_circle_sizes, pch = 16, col = "black")  # Add circles on top of base map

# Step 8: Add a custom legend to show the quarter_numeric values and corresponding circle sizes
legend("topright", legend = unique_quarters, 
       pch = 16, 
       pt.cex = circle_sizes,  # Show unique circle sizes in the legend
       col = "black", 
       title = "Quarters (Circle Size)")

Sagaing_QuarterYear_PPP_Numeric <- ppp(
  x = coords[, 1],  # Longitude (x-coordinates)
  y = coords[, 2],  # Latitude (y-coordinates)
  window = Sagaing_Owin,  # Spatial window (Sagaing_Owin)
  marks = Sagaing_QuarterYear_Sf$quarter_numeric  # Use numeric for KDE analysis (time-based)
)

Sagaing_QuarterYear_KDE <- spattemp.density(Sagaing_QuarterYear_PPP_Numeric)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
multi_color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))


# Define the quarters for which you want to plot the KDE
tims <- c(20211, 20212, 20213, 20214, 20221, 20222, 20223, 20224, 
          20231, 20232, 20233, 20234, 20241, 20242)

#Set up the plotting window to display multiple plots in a grid (4 rows, 4 columns)
par(mfrow=c(4, 4), mar=c(2, 2, 2, 2), oma=c(0, 0, 2, 0))  # Adjust margins and add outer margin for title

# Loop through the 'tims' vector and plot the KDE for each time point (from left to right)
for(i in tims) { 
  plot(Sagaing_QuarterYear_KDE, i, 
       override.par=FALSE,  # Keep the graphical parameters unchanged
       fix.range=TRUE,      # Fix the range of the KDE
       main=paste("Quarter", i),  # Title for each plot
       ribside = "right", 
       col = multi_color_palette(100))  # Use 'inferno' color gradient
  
  # Remove the x and y axis labels for a cleaner look
  axis(1, labels=FALSE)
  axis(2, labels=FALSE)
}

#dd a unified title across all plots (e.g., the duration of analysis)
mtext("Quarterly Spatio-Temporal KDE for Civilian Conflict Intensity in Sagaing", 
      outer=TRUE, cex=1.5, font=2)

multi_color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))

# Save the animation as a GIF
saveGIF({
  # Loop over the valid quarter times
  for(i in tims){ 
    suppressWarnings({
      plot(Sagaing_QuarterYear_KDE, i, 
           override.par=FALSE,  # Keep graphical parameters unchanged
           fix.range=TRUE,      # Fix the range of the KDE
           main=paste("Sagaing Civilian Related Conflict KDE at Quarter", i),  # Title for each plot
           ribside = "right",  # Legend on the right
           col = multi_color_palette(100))  # Apply multi-color gradient
    })
  }
}, movie.name = "Sagaingkde_animation.gif", interval = 0.5, ani.width = 800, ani.height = 800)

Spatio-Temporal KDE analysis for Sagaing across multiple quarters, we observe significant spatial clustering of civilian-related conflicts over time. The conflict intensity remains localized in specific regions, with fluctuations in density between quarters.

Key Observations:

  • Consistent Clustering: Conflict events are consistently clustered in the central and southern regions of Sagaing, as indicated by the green to red areas over time. These regions show recurring intensity, highlighting persistent conflict zones.

  • Temporal Fluctuations: Some quarters, such as 20212, 20213, and 20222, exhibit relatively lower conflict intensity, while 20224 shows a resurgence in intensity.

  • Localized Hotspots: The southern region shows a spike in density during 20241, followed by a steady reduction in the later quarters.

Conclusion:

There is ongoing clustering in the same geographic regions of Sagaing, with peaks of conflict intensity in specific quarters. This indicates that while the overall distribution of conflicts remains stable, certain periods witness heightened violence.

10.0 Perform 2nd-Order Spatio-Temporal Point Patterns Analysis:

In this analysis, we aim to determine whether civilian-related conflict events in States are clustered or randomly distributed across both space and time. To achieve this, we use spatial point pattern analysis, specifically the K-function, to explore clustering between conflict events over different years. This approach helps detect patterns of clustering or dispersion across varying spatial scales, comparing conflict events year by year.

The following steps were taken:

  1. Creation of the Point Pattern Object:

    • The spatial point pattern object (Sagaing_ppp_years) was created using the conflict event data for Sagaing, where the coordinates (longitude and latitude) were included along with the quarter-year information as marks to differentiate between events from different years.

  2. Rescaling:

    • The point pattern object was rescaled to kilometers to ensure that the distance measurements in the analysis were meaningful and consistent with spatial interpretation.

  3. Cross-K Function Computation:

    • For each pair of years (2021 vs. 2022, 2021 vs. 2023, 2022 vs. 2023, etc.), the cross-K function was computed. This function measures how events in one year are spatially related to events in another year. The results from each comparison were stored in the {state}_k_results list for further interpretation.

10.1 Yangon

For the conflict events in Sagaing, we performed a detailed spatial point pattern analysis by computing the cross-K function for each pair of years (e.g., 2021 vs. 2022, 2021 vs. 2023, etc.). The cross-K function allows us to investigate how conflict events in one year relate to conflict events in another year, revealing whether the events are clustered, dispersed, or randomly distributed across space between different time periods.

Null and Alternative Hypotheses:

  • Ho (Null Hypothesis): The distribution of conflict events in Yangon follows Complete Spatio-Temporal Randomness (CSTR), meaning the events are randomly distributed in both space and time.

  • H1 (Alternative Hypothesis): The distribution of conflict events in Yangon does not follow CSTR and instead shows signs of clustering or dispersion in both space and time.

10.1.1 Yangon’s Cross K Calculation

# Create the spatial point pattern (ppp) object for Yangon
Yangon_ppp_years <- ppp(
  x = Yangon_ACLED_co[, 1],  # Longitude
  y = Yangon_ACLED_co[, 2],  # Latitude
  window = Yangon_Owin  # Spatial observation window (owin object)
)

# Add a mark that groups quarters by year (e.g., 2021, 2022, 2023, 2024)
Yangon_ppp_years$marks <- as.factor(substr(Yangon_QuarterYear_Sf$quarter_numeric, 1, 4))
Yangon_ppp_years <- rescale(Yangon_ppp_years, 1000, "km")
# Store the unique years in the data
unique_years <- levels(Yangon_ppp_years$marks)

# Function to compute cross-K function between two years
cross_k_function_by_year <- function(ppp_object, year1, year2) {
  # Compute the cross-K function for the two selected years within the same ppp object
  K_cross <- Kcross(ppp_object, i = year1, j = year2)
  return(K_cross)
}

# Initialize an empty list to store each year's comparison
Yangon_k_results <- list()

# Loop through pairs of years to compute and store cross-K functions
for (i in 1:(length(unique_years) - 1)) {
  for (j in (i + 1):length(unique_years)) {
    year1 <- unique_years[i]
    year2 <- unique_years[j]
    
    # Print progress
    cat("Computing K-Function for", year1, "and", year2, "\n")
    
    # Compute the cross-K function for these two years and store it in a named variable
    result_name <- paste0("Yangon_k_", year1, "_vs_", year2)
    Yangon_k_results[[result_name]] <- cross_k_function_by_year(Yangon_ppp_years, year1, year2)
  }
}

# Save individual year comparison results as variables
Yangon_k_2021_vs_2022 <- Yangon_k_results[["Yangon_k_2021_vs_2022"]]
Yangon_k_2021_vs_2023 <- Yangon_k_results[["Yangon_k_2021_vs_2023"]]
Yangon_k_2021_vs_2024 <- Yangon_k_results[["Yangon_k_2021_vs_2024"]]
Yangon_k_2022_vs_2023 <- Yangon_k_results[["Yangon_k_2022_vs_2023"]]
Yangon_k_2022_vs_2024 <- Yangon_k_results[["Yangon_k_2022_vs_2024"]]
Yangon_k_2023_vs_2024 <- Yangon_k_results[["Yangon_k_2023_vs_2024"]]

# Set up the layout for multiple plots (2 rows, 3 columns for all year comparisons)
par(mfrow = c(2, 1))
# Plotting each comparison
plot(Yangon_k_2021_vs_2022, main = "K-Function: Yangon Civilian Related Conflict KDE 2021 vs 2022")
plot(Yangon_k_2021_vs_2023, main = "K-Function: Yangon Civilian Related Conflict KDE 2021 vs 2023")

plot(Yangon_k_2021_vs_2024, main = "K-Function: Yangon Civilian Related Conflict KDE 2021 vs 2024")
plot(Yangon_k_2022_vs_2023, main = "K-Function: Yangon Civilian Related Conflict KDE 2022 vs 2023")

plot(Yangon_k_2022_vs_2024, main = "K-Function: Yangon Civilian Related Conflict KDE 2022 vs 2024")
plot(Yangon_k_2023_vs_2024, main = "K-Function: Yangon Civilian Related Conflict KDE 2023 vs 2024")

10.1.2 Yangon’s Cross K Analysis

10.1.2.1 Yangon’s Individual Years Comparision

10.1.2.1.1 2021 vs 2022

  • The observed K-function lies consistently above the CSR envelope, indicating significant clustering between events from 2021 and 2022 at distances up to approximately 25 km.

  • The clustering effect becomes more prominent at larger distances, suggesting that events from 2021 tend to cluster with events from 2022 across larger spatial scales.

10.1.2.1.2 2021 vs 2023

  • The K-function for 2021 vs 2023 is also above the CSR envelope, indicating clustering between these two years.

  • However, the clustering is slightly less prominent compared to 2021 vs 2022, suggesting that while clustering persists, the relationship between events from 2021 and 2023 is somewhat weaker than between 2021 and 2022.

10.1.2.1.3 2021 vs 2024

  • The clustering is still evident between 2021 and 2024, but the clustering intensity appears to be weaker at larger distances. The observed K-function shows less deviation from the CSR envelope, indicating that the spatial clustering effect is declining as we compare events separated by more time.
10.1.2.1.4 2022 vs 2023

  • The K-function here indicates significant clustering at smaller distances, but the intensity decreases more rapidly beyond 15 km.

  • The spatial clustering between 2022 and 2023 is still present but somewhat weaker compared to other year-pair comparisons.

10.1.2.1.5 2022 vs 2024

  • The 2022 vs 2024 K-function shows clustering, though it is less pronounced than in the previous year comparisons.

  • The clustering effect diminishes as the events from these two years become more spatially dispersed, particularly beyond 15 km.

10.1.2.1.6 2023 vs 2024

  • The K-function for 2023 vs 2024 shows the least clustering effect of all the comparisons.

  • There is minimal deviation from the CSR line, indicating that the events from these two years are closer to being randomly distributed, with only mild clustering observed at smaller distances.

10.1.2.1 Overall Years Comparison

  • Clustering is strongest between events that are closer in time (e.g., 2021 vs 2022), with the strength of clustering generally decreasing as the time difference between events increases (e.g., 2021 vs 2024).

  • This suggests that civilian conflict events in Yangon have a spatial-temporal dependence, where events in one year are more likely to occur near events from the following year, but this effect diminishes over time.

10.2 Mandalay

For the conflict events in Mandalay, we performed a detailed spatial point pattern analysis by computing the cross-K function for each pair of years (e.g., 2021 vs. 2022, 2021 vs. 2023, etc.). The cross-K function allows us to investigate how conflict events in one year relate to conflict events in another year, revealing whether the events are clustered, dispersed, or randomly distributed across space between different time periods.

Null and Alternative Hypotheses:

  • Ho (Null Hypothesis): The distribution of conflict events in Mandalay follows Complete Spatio-Temporal Randomness (CSTR), meaning the events are randomly distributed in both space and time.

  • H1 (Alternative Hypothesis): The distribution of conflict events in Mandalay does not follow CSTR and instead shows signs of clustering or dispersion in both space and time.

10.2.1 Mandalay’s Cross K Calculation

# Create the spatial point pattern (ppp) object for Mandalay
Mandalay_ppp_years <- ppp(
  x = Mandalay_ACLED_co[, 1],  # Longitude
  y = Mandalay_ACLED_co[, 2],  # Latitude
  window = Mandalay_Owin  # Spatial observation window (owin object)
)

# Add a mark that groups quarters by year (e.g., 2021, 2022, 2023, 2024)
Mandalay_ppp_years$marks <- as.factor(substr(Mandalay_QuarterYear_Sf$quarter_numeric, 1, 4))
# Store the unique years in the data
unique_years <- levels(Mandalay_ppp_years$marks)
Mandalay_ppp_years <- rescale(Mandalay_ppp_years, 1000, "km")

# Function to compute cross-K function between two years
cross_k_function_by_year <- function(ppp_object, year1, year2) {
  # Compute the cross-K function for the two selected years within the same ppp object
  K_cross <- Kcross(ppp_object, i = year1, j = year2)
  return(K_cross)
}

# Initialize an empty list to store each year's comparison
Mandalay_k_results <- list()

# Loop through pairs of years to compute and store cross-K functions
for (i in 1:(length(unique_years) - 1)) {
  for (j in (i + 1):length(unique_years)) {
    year1 <- unique_years[i]
    year2 <- unique_years[j]
    
    # Print progress
    cat("Computing K-Function for", year1, "and", year2, "\n")
    
    # Compute the cross-K function for these two years and store it in a named variable
    result_name <- paste0("Mandalay_k_", year1, "_vs_", year2)
    Mandalay_k_results[[result_name]] <- cross_k_function_by_year(Mandalay_ppp_years, year1, year2)
  }
}

# Save individual year comparison results as variables
Mandalay_k_2021_vs_2022 <- Mandalay_k_results[["Mandalay_k_2021_vs_2022"]]
Mandalay_k_2021_vs_2023 <- Mandalay_k_results[["Mandalay_k_2021_vs_2023"]]
Mandalay_k_2021_vs_2024 <- Mandalay_k_results[["Mandalay_k_2021_vs_2024"]]
Mandalay_k_2022_vs_2023 <- Mandalay_k_results[["Mandalay_k_2022_vs_2023"]]
Mandalay_k_2022_vs_2024 <- Mandalay_k_results[["Mandalay_k_2022_vs_2024"]]
Mandalay_k_2023_vs_2024 <- Mandalay_k_results[["Mandalay_k_2023_vs_2024"]]
# Set up the layout for multiple plots (2 rows, 3 columns for all year comparisons)
par(mfrow = c(2, 1))

# Plotting each comparison
plot(Sagaing_k_2021_vs_2022, main = "K-Function: Sagaing Civilian Related Conflict KDE 2021 vs 2022")
plot(Sagaing_k_2021_vs_2023, main = "K-Function: Sagaing Civilian Related Conflict KDE 2021 vs 2023")

plot(Sagaing_k_2021_vs_2024, main = "K-Function: Sagaing Civilian Related Conflict KDE 2021 vs 2024")
plot(Sagaing_k_2022_vs_2023, main = "K-Function: Sagaing Civilian Related Conflict KDE 2022 vs 2023")

plot(Sagaing_k_2022_vs_2024, main = "K-Function: Sagaing Civilian Related Conflict KDE 2022 vs 2024")
plot(Sagaing_k_2023_vs_2024, main = "K-Function: Sagaing Civilian Related Conflict KDE 2023 vs 2024")

10.2.2 Mandalay’s Cross K Analysis

10.2.2.1 Mandalay’s Individual Years Comparison

10.1.2.1.1 2021 vs 2022

  • The K-function for 2021 vs 2022 shows significant clustering at all distances. The observed K-function is consistently above the CSR envelope, indicating that conflict events from 2021 tend to cluster near those from 2022.
  • The clustering effect appears to intensify as the distance increases, suggesting a broad spatial clustering pattern for conflicts in 2021 and 2022.
10.2.2.1.2 2021 vs 2023

  • the K-function for 2021 vs 2023 also shows clustering, but it is slightly less pronounced than in the 2021 vs 2022 comparison.
  • The clustering becomes more noticeable at distances beyond 10 km, implying that conflict events from these two years tend to cluster over medium to larger distances.
  • However, the clustering effect weakens somewhat compared to the previous year comparison.
10.2.2.1.3 2021 vs 2024

  • Clustering between 2021 and 2024 is still present, but the intensity is weaker compared to the previous year pairs.
  • The observed K-function does show some clustering above the CSR line, but the deviation from randomness is smaller, suggesting that events from 2021 and 2024 are less spatially related than the closer years (2021 vs 2022).
10.2.2.1.4 2022 vs 2023

  • The clustering between 2022 and 2023 remains significant, especially at larger distances, with the observed K-function clearly above the CSR envelope.
  • Similar to earlier comparisons, this suggests that conflict events from these two years cluster over wider distances, continuing the spatial clustering trend seen in other comparisons.
10.2.2.5 2022 vs 2024

  • The 2022 vs 2024 K-function shows clustering, although the intensity is weaker, particularly at smaller distances.
  • There is a slight increase in clustering at larger distances, but the deviation from the CSR envelope is less pronounced compared to the earlier year comparisons.
10.2.2.6 2023 vs 2024

  • The K-function for 2023 vs 2024 shows clustering, though this is the weakest of all the comparisons.

  • The clustering effect becomes more noticeable at distances greater than 10 km, but overall, the spatial relationship between events from 2023 and 2024 is closer to random compared to earlier year comparisons.

10.2.2.1 Mandalay Overall Years Comparison

  • Similar to the Yangon analysis, the spatial clustering is strongest between conflict events that are closer in time, such as 2021 vs 2022 and 2022 vs 2023.
  • As the time gap between years increases, the clustering effect becomes weaker, with the least clustering observed between 2023 and 2024.

These results suggest that conflicts in Mandalay show spatio-temporal clustering, with events from one year often clustering near those from the next, but this effect diminishes as the time gap between years grows.

10.3 Sagaing

For the conflict events in Sagaing, we performed a detailed spatial point pattern analysis by computing the cross-K function for each pair of years (e.g., 2021 vs. 2022, 2021 vs. 2023, etc.). The cross-K function allows us to investigate how conflict events in one year relate to conflict events in another year, revealing whether the events are clustered, dispersed, or randomly distributed across space between different time periods.

Null and Alternative Hypotheses:

  • Ho (Null Hypothesis): The distribution of conflict events in Sagaing follows Complete Spatio-Temporal Randomness (CSTR), meaning the events are randomly distributed in both space and time.

  • H1 (Alternative Hypothesis): The distribution of conflict events in Sagaing does not follow CSTR and instead shows signs of clustering or dispersion in both space and time.

10.3.1 Sagaing’s Cross K Calculation

Sagaing_ppp_years <- ppp(
  x = Sagaing_ACLED_co[, 1],  # Longitude
  y = Sagaing_ACLED_co[, 2],  # Latitude
  window = Sagaing_Owin  # Spatial observation window (owin object)
)

# Add a mark that groups quarters by year (e.g., 2021, 2022, 2023, 2024)
Sagaing_ppp_years$marks <- as.factor(substr(Sagaing_QuarterYear_Sf$quarter_numeric, 1, 4))
Sagaing_ppp_years <- rescale(Sagaing_ppp_years, 1000, "km")
# Store the unique years in the data
unique_years <- levels(Sagaing_ppp_years$marks)

# Function to compute cross-K function between two years
cross_k_function_by_year <- function(ppp_object, year1, year2) {
  # Compute the cross-K function for the two selected years within the same ppp object
  K_cross <- Kcross(ppp_object, i = year1, j = year2)
  return(K_cross)
}

# Initialize an empty list to store each year's comparison
Sagaing_k_results <- list()

# Loop through pairs of years to compute and store cross-K functions
for (i in 1:(length(unique_years) - 1)) {
  for (j in (i + 1):length(unique_years)) {
    year1 <- unique_years[i]
    year2 <- unique_years[j]
    
    # Print progress
    cat("Computing K-Function for", year1, "and", year2, "\n")
    
    # Compute the cross-K function for these two years and store it in a named variable
    result_name <- paste0("Sagaing_k_", year1, "_vs_", year2)
    Sagaing_k_results[[result_name]] <- cross_k_function_by_year(Sagaing_ppp_years, year1, year2)
  }
}

# Save individual year comparison results as variables
Sagaing_k_2021_vs_2022 <- Sagaing_k_results[["Sagaing_k_2021_vs_2022"]]
Sagaing_k_2021_vs_2023 <- Sagaing_k_results[["Sagaing_k_2021_vs_2023"]]
Sagaing_k_2021_vs_2024 <- Sagaing_k_results[["Sagaing_k_2021_vs_2024"]]
Sagaing_k_2022_vs_2023 <- Sagaing_k_results[["Sagaing_k_2022_vs_2023"]]
Sagaing_k_2022_vs_2024 <- Sagaing_k_results[["Sagaing_k_2022_vs_2024"]]
Sagaing_k_2023_vs_2024 <- Sagaing_k_results[["Sagaing_k_2023_vs_2024"]]

# Save individual year comparison results as variables
Sagaing_k_2021_vs_2022 <- Sagaing_k_results[["Sagaing_k_2021_vs_2022"]]
Sagaing_k_2021_vs_2023 <- Sagaing_k_results[["Sagaing_k_2021_vs_2023"]]
Sagaing_k_2021_vs_2024 <- Sagaing_k_results[["Sagaing_k_2021_vs_2024"]]
Sagaing_k_2022_vs_2023 <- Sagaing_k_results[["Sagaing_k_2022_vs_2023"]]
Sagaing_k_2022_vs_2024 <- Sagaing_k_results[["Sagaing_k_2022_vs_2024"]]
Sagaing_k_2023_vs_2024 <- Sagaing_k_results[["Sagaing_k_2023_vs_2024"]]
par(mfrow = c(2, 1))
# Plotting each comparison
plot(Sagaing_k_2021_vs_2022, main = "K-Function: Sagaing Civilian Related Conflict KDE 2021 vs 2022")
plot(Sagaing_k_2021_vs_2023, main = "K-Function: Sagaing Civilian Related Conflict KDE 2021 vs 2023")

plot(Sagaing_k_2021_vs_2024, main = "K-Function: Sagaing Civilian Related Conflict KDE 2021 vs 2024")
plot(Sagaing_k_2022_vs_2023, main = "K-Function: Sagaing Civilian Related Conflict KDE 2022 vs 2023")

plot(Sagaing_k_2022_vs_2024, main = "K-Function: Sagaing Civilian Related Conflict KDE 2022 vs 2024")
plot(Sagaing_k_2023_vs_2024, main = "K-Function: Sagaing Civilian Related Conflict KDE 2023 vs 2024")

10.3.2 Sagaing’s Cross K Analysis

10.3.2.1 Sagaing’s Individual Years Comparison

10.3.2.1.1 2021 vs 2022

  • The observed K-function r is consistently higher than the CSR (Complete Spatial Randomness) envelope across almost all distances. This indicates that the conflict events in 2021 and 2022 are significantly clustered in space. The clustering effect becomes more pronounced at larger distances, especially after 20 km.

  • The trans K-function also follows a similar pattern, showing clustering beyond the expected CSR pattern.

  • The Kbord function demonstrates that boundary corrections had minimal impact on the general clustering trends observed.

10.3.2.1.2 2021 vs 2023

  • Similar to the comparison between 2021 and 2022, the K-function in 2021 vs 2023 also shows a high degree of clustering, with obs(r) being higher than CSR across nearly all distances.

  • There is a slightly more pronounced clustering effect around 30 to 40 km, where the gap between the observed K and CSR becomes wider. This suggests an increase in conflict density over time.

10.3.2.1.3 2021 vs 2024

  • the K-function for 2021 vs 2024 shows a consistent pattern of clustering across distances, although the observed values tend to rise more sharply after 40 km, indicating that conflicts in 2024 were more dispersed initially but clustered significantly at larger spatial scales.

  • The results suggest that the spatial distribution of conflicts in 2024 is slightly different, possibly influenced by external factors that caused conflicts to be more spread out in space.

10.3.2.1.4 2022 vs 2023

  • The comparison between 2022 and 2024 reveals strong clustering over a wide range of distances, with significant increases after 30 km.
  • This suggests that conflicts in 2024 are more spatially concentrated at larger distances compared to 2022.
10.3.2.1.5 2022 vs 2024

  • The comparison between 2022 and 2024 reveals strong clustering over a wide range of distances, with significant increases after 30 km.

  • This suggests that conflicts in 2024 are more spatially concentrated at larger distances compared to 2022.

10.3.2.1.6 2023 vs 2024

  • The clustering trend continues between 2023 and 2024, with notable increases in clustering beyond 40 km. The gap between the observed K and CSR remains large, indicating that the conflict events in 2023 and 2024 were not randomly distributed but exhibited significant spatial clustering.

10.3.2.1 Sagaing Overall Years Comparison

  • Across all years (2021-2024), the K-functions suggest that civilian-related conflict events in Sagaing are significantly clustered rather than randomly distributed. The degree of clustering varies slightly across the years, with some years (such as 2024) showing stronger clustering at larger distances (over 40 km). This indicates that the conflict dynamics in Sagaing are evolving over time, possibly due to changes in political, social, or environmental factors.

11.0 Conclusion

This spatio-temporal analysis of civilian-related conflicts in Yangon, Sagaing, and Mandalay from 2021 to 2024 shows consistent clustering of conflict events, suggesting that conflicts are not randomly distributed but rather concentrated in specific hotspots across all three states. The K-function comparisons between years indicate that conflicts in these regions are becoming increasingly clustered, especially as distances between events grow larger. This analysis highlights the persistent and expanding nature of conflict zones in Myanmar.

However, this analysis only focuses on the number of conflicts and fatalities within civilian-related conflict events. A more comprehensive understanding could be achieved by introducing categorical analysis, examining each conflict’s specific causes, involved groups, and other characteristics to uncover deeper insights into the dynamics driving these conflict patterns.

Areas to Avoid:

  1. Yangon: Conflict events in Yangon are significantly clustered, particularly in the city center and the outskirts towards the east and south. These areas have shown increasing conflict over time, with 2023 and 2024 seeing higher concentrations of incidents.

  2. Sagaing: The northern regions of Sagaing, particularly towards the Sagaing Hills and nearby rural settlements, exhibit high clustering of conflicts. Avoid traveling near these conflict hotspots, especially as the spatial extent of clustering has grown between 2022 and 2024.

  3. Mandalay: In Mandalay, the central urban areas and southern outskirts demonstrate increased clustering of conflict events. The comparison between 2021 and 2024 reveals that conflict events are more widespread in these areas, indicating a need for caution in the central and southern parts of the state.

By addressing both spatial and temporal patterns, authorities and humanitarian organizations can better strategize responses and interventions in these conflict-prone areas.

12.0 References

13.0 Miscellaneous for Usability

In this section, we ensure that all the critical data, variables, and results generated throughout the analysis are exported and saved for subsequent usage. By saving everything in .rds format, we create an easy-to-load, reusable dataset and environment that allows future analysis without the need to re-run the entire process. This step enhances reproducibility and usability for further research or adjustments.

save_dir <- "rds_files/"
if (!dir.exists(save_dir)) {
  dir.create(save_dir)
}

vars <- ls()

# Save each variable as an RDS file
for (var in vars) {
  saveRDS(get(var), file = paste0(save_dir, var, ".rds"))
}